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INTRODUCTION 


Models  used  to  simulate  the  inherent  capabilities  of  modem  defense  systems  are 
employed  extensively  by  today's  analyst.  Because  of  rising  costs  on  actual  hardware 
testing  and  training,  detailed  models  that  provide  reliable  output  are  routinely  sought  after 
or  developed.  This  is  particularly  true  in  the  area  of  survivability  analysis. 

Combat  survivability,  by  definition,  is  the  ability  of  a  vehicle  to  avoid  and/or 
withstand  a  man-made  hostile  environment.  The  threat  in  such  an  environment  consists  of 
radars,  guns,  missiles,  and  so  on.  In  survivability  analysis,  the  effects  of  each  threat  in 
this  hostile  environment  can  be  studied  individually  using  sophisticated  models.  In  aircraft 
combat  survivability,  a  necessary  input  to  these  models  is  a  position  and  attitude  history  of 
the  aerodynamic  vehicle  being  analyzed.  The  generation  of  such  a  flight  path  for  the  vehicle 
is  therefore  a  necessary  first  step  toward  assessing  the  susceptibility  and/or  the  vulnerability 
of  the  vehicle  to  a  particular  threat  system. 

MOFO  (Model  Of  Flying  Objects)  is  a  computer  model  that  generates  flight  paths  of 
aerodynamic  vehicles.  It  can  simulate  conventional  aircraft  and,  with  minor  modifications, 
bank-to-tum  missiles.  The  user  provides  both  the  necessary  aerodynamic  data  and  power 
parameters  of  the  vehicle  in  one  file  and  the  commands  that  specify  the  intended  flight  path 
in  another  file.  The  model  computes  a  six-degree-of-freedom  trajectory  consisting  of  a  time 
history  of  position  and  attitude,  using  equations  found  in  [1].  Such  output  can  be  used  as 
input  to  models  of  radar  detection  or  missile  flyout  as  described  in  [2], 

The  purpose  of  this  manual  is  to  provide  the  user  with  the  critical  information 
necessary  to  set  up  and  run  the  model.  The  first  section  provides  a  functional  description 
of  MOFO.  It  delineates  the  computational  requirements  and  the  model’s  coordinate 
systems.  A  second  section  covers  setup  and  run  instructions,  and  includes  sample  files  and 
their  respective  formats.  The  manual  concludes  with  improvements  that  should  be 
implemented  in  the  model  at  a  later  date. 
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FUNCTIONAL  DESCRIPTION 


COMPUTATIONAL  DETAILS 

MOFO  is  written  in  FORTRAN  77  and  currently  runs  on  a  Micro  VAX  II  computer. 
The  present  model  requires  roughly  350  blocks  of  memory  on  a  Digital  RA60  disk  drive. 
This  number  includes  all  source  code  files,  text  and  object  libraries,  the  executable  image, 
and  ".COM”  files.  It  does  not  include  the  storage  requirement  for  either  the  input  or  output 
files,  which  is  a  function  of  the  user's  analytical  objective  and  disk  capacity.  Except  for  a 
few  long  variable  names  that  are  allowed  in  the  VAX/VMS  FORTRAN  language,  standard 
coding  conventions  have  been  used.  The  code  itself  is  unclassified,  but  the  user  may  have 
classified  input  parameters  to  use.  These  attributes  should  serve  to  increase  the 
transportability  of  MOFO  to  other  machine  environments. 

In  its  present  form,  MOFO  contains  1,300  lines  of  code  divided  into  23  separate 
modules:  subroutines,  functions,  and  the  main  program.  Figure  1  provides  an  overview 
of  the  model's  structure  using  a  call  tree  diagram.  Implicit  double  precision  is  used 
throughout  in  an  effort  to  minimize  round-off  errors.  Common  blocks  are  placed  in  a 
separate  file  and  an  "INCLUDE"  statement  is  used  by  the  appropriate  modules  to 
incorporate  the  variable  transfer  structure.  More  specific  information  contained  in 
appendix  A  describes  both  variable  definitions  and  details  of  each  module's  function. 
MOFO  is  internally  documented  as  well;  the  source  code  is  contained  in  appendix  B. 

Whether  the  user  is  simply  modifying  a  module  of  the  model  or  creating  the 
executable  image  from  scratch,  the  CPU  time  required  to  compile  and  link  the  program  is 
nominal.  The  user  should  structure  the  MOFO  model  directory  as  described  in  the  file 
MODMOFO.COM  (a  hardcopy  is  located  in  appendix  C).  With  this  directory  structure 
intact,  the  MODMOFO.COM  file  may  be  used  to  minimize  the  effort  involved  in  the 
modification  process.  Actual  run  times  for  the  model  are,  of  course,  dependent  on  the 
length  of  the  scenario  for  which  the  model  is  trying  to  generate  a  trajectory.  But  the  ratio  of 
CPU  time  to  real  time  is  approximately  one  to  one. 

INERTIAL  REFERENCE  SYSTEM 

MOFO  models  the  flight  of  a  single  aerodynamic  vehicle  as  it  attempts  to  follow  a 
commanded  route  via  waypoints.  Waypoints  are  positions  in  three-dimensional  space  that 
the  aircraft  attempts  to  fly  through  in  an  effort  to  produce  the  desired  trajectory.  The 
trajectory  is  generated  in  an  inertial  reference  frame  fixed  on  the  surface  of  a  flat  earth. 
Position  is  described  via  Cartesian  coordinates  X,  Y,  and  Z  in  this  frame.  These 
coordinates  are  measured  in  meters.  As  shown  in  figure  2,  the  positive  X-axis  is  directed 
due  east,  the  positive  Y-axis  due  north,  and  the  positive  Z-axis  upwards.  The  earth's 
atmosphere  in  MOFO  is  assumed  stationary  and  homogeneous.  The  atmospheric  pressure, 
density,  and  temperature  profiles  are  those  of  the  U.S.  standard  atmosphere.  These 
profiles  are  computed  from  expressions  found  in  [3]. 

LOCAL  HORIZON 

While  MOFO  is  generating  a  flight  path  of  positional  data,  it  is  also  keeping  track  of 
the  vehicle's  attitude  relative  to  the  local  horizon  reference  frame.  The  local  horizon  is  a 
right-handed  system  with  its  origin  based  at  the  aircraft's  center  of  gravity  (c.g.).  Each 
axis  is  parallel  to  its  corresponding  axis  in  the  inertial  reference  frame  (figure  2).  To 
specify  the  aerodynamic  forces  and  the  attitude  of  the  aircraft,  two  additional  reference 
frames  are  needed:  the  wind  system  and  the  body  system. 
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Figure  1.  Call  tree  for  MOFO  modules 


Figure  2.  MOFO  external  coordinate  system 

WIND  SYSTEM 

The  wind  system  specifies  the  direction  of  the  aircraft's  velocity  vector  and  therefore 
serves  as  a  natural  frame  of  reference  for  describing  the  aerodynamic  forces.  Three 
rotations  are  needed  to  arrive  at  the  wind  system  from  the  local  horizon  (see  figure  3). 
Consider  an  aircraft  that  is  flying  straight  and  level  with  its  velocity  vector  pointing 
northward  along  the  yh-axis  and  its  starboard  wing  pointing  eastward  along  the  Xh-axis. 
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Figure  3.  Rotations  from  the  local  horizon  to  the  wind  system 


The  following  sequence  of  rotations  about  the  c.g.  is  performed: 

*  The  velocity  vector  is  rotated  in  a  left-handed  sense  through  an  angle  % 
about  the  zh-axis  in  the  xh  yh  plane  to  arrive  at  the  intermediate  xi  yi  z\ 
system.  The  angle  %  is  the  velocity  heading  relative  to  north. 

*  The  velocity  vector  is  then  rotated  in  a  right-handed  sense  through  the  angle 
7  about  the  xi-axis  in  the  yi  z\  plane  to  arrive  at  the  intermediate  X2  y2  Z2 
system.  The  angle  7  is  the  velocity  pitch  angle. 

*  Last,  the  aircraft  is  rotated  in  a  right-handed  sense  through  angle  4  about 
the  y2-axis  (velocity  vector)  in  the  X2  Z2  plane  to  arrive  at  the  xw  yw  zw  wind 
system.  The  angle  4  is  the  velocity  roll  angle. 
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BODY  SYSTEM 

The  body  system  is  used  to  specify  the  orientation  of  the  aircraft  relative  to  the 
oncoming  wind.  Two  angular  rotations  are  needed  to  arrive  at  the  body  system  from  the 
wind  system  (see  figure  4).  The  following  sequence  of  rotations  about  the  c.g.  is 
performed: 


•  The  aircraft  is  rotated  in  a  left-handed  sense  through  an  angle  ct  about  the 
zw-axis  in  the  xw  yw  plane  to  arrive  at  the  intermediate  X3  y3  Z3  system.  The 
angle  a  is  the  sideslip  angle. 

•  The  aircraft  is  then  rotated  in  a  right-handed  sense  through  an  angle  a 
about  the  x3-axis  in  the  y3  Z3  plane  to  arrive  at  the  xb  yb  zb  body  system. 
The  angle  a  is  the  angle  of  attack. 


Figure  4.  Rotations  from  the  wind  system  to  the  body  system 


Because  the  aircraft  is  assumed  to  fly  without  sideslip,  o  will  always  equal  zero  in 
MOFO.  All  the  other  angles  (%,  y,  ji,  and  a)  are  written  to  an  output  file.  These  angles 
appear  in  the  equations  of  motion  as  found  in  [1], 
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The  body  system  can  also  be  achieved  directly  from  the  local  horizon  in  a  manner 
analogous  to  the  transformation  to  the  wind  system  from  the  local  horizon.  The  velocity 
angles  in  figure  3  (%,  y,  and  |i)  are  replaced  by  the  body  angles  (Xb,  Yb,  and  jib, 
respectively).  These  angles  are  written  to  an  output  file  along  with  the  aircraft's  positional 
data.  This  file  is  separate  from  the  file  that  contains  the  velocity  angles  and  the  angle  of 
attack.  It  is  this  file  that  would  be  used  as  input  to  other  models  requiring  six-degree-of- 
ffeedom  information.  For  example,  a  detailed  radar  model  would  need  position  and  attitude 
information  to  determine  the  appropriate  radar  cross  section  (RCS)  value  along  the  radar's 
line  of  sight. 

EXTERNAL  VERSUS  INTERNAL  COORDINATE  SYSTEMS 

The  external  coordinate  systems  just  described  are  used  to  ensure  the  compatibility  of 
MOFO  with  other  models  used  in  the  Surface  and  Air  Combat  Systems  Program.  The 
coordinate  systems  used  internally  in  MOFO  are  analogous  but  slightly  different  (see  [1]). 
For  instance,  the  following  axis  inversions  are  required  to  transform  the  external  inertial 
system  to  the  internal  inertial  system: 


Y,  =  -Ye  (1) 

and 

Z/  =  -Zr  (2) 


where  /  denotes  the  internal  system,  and  E  denotes  the  external  system  as  described 
under  "Inertial  Reference  System."  The  subscript  E  was  dropped  before  for  simplicity  of 
notation.  These  axis  inversions  create  a  right-handed  coordinate  system  in  the  internal 
system.  Because  internal  headings  are  measured  from  the  inertial  Xj-axis  and  external 
headings  are  from  the  inertial  Yg-axis,  a  third  transformation  is: 


X/  =  %E  -  90°  .  (3) 

Before  the  positional  and  angular  data  are  written  to  the  output  files,  the  appropriate 
variables  are  transformed  back  to  the  external  coordinate  system  for  consistency.  The  user 
should  be  aware  of  these  differences  between  the  two  systems  in  the  event  that  later 
modifications  are  required. 
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SETUP  AND  RUN  INSTRUCTIONS 


DCL  COMMAND  RUN  FILE 

On  the  assumption  that  the  user  is  familiar  with  running  a  ".COM"  file  interactively 
or  in  batch  mode,  the  command  file  in  the  VAX/VMS  operating  system  for  running  the 
MOFO  model  should  resemble: 

$  ASSIGN  AIRCRAFT. IN  FOROOl 
$  ASSIGN  COMMAND. IN  FOR002 
$  ! 

$  ASSIGN  POSITION. OUT  FOR007 
$  ASSIGN  VELOCITY. OUT  FOR008 
$  ASSIGN  ERROR. OUT  FOR009 
$  ! 

$  RUN  MOFO.EXE 
$  ! 

$  DEASSIGN  FOROOl 
$  DEASSIGN  FOR002 
$  DEASSIGN  FOR007 
$  DEASSIGN  FOR008 
$  DEASSIGN  FOR009 
$  ! 

$  EXIT 

Logical  units  1  and  2  are  the  aircraft  characteristics  file  and  command  file, 
respectively.  Logical  unit  9  contains  error  messages  that  alert  the  user  to  possible 
problems.  Logical  unit  7  is  the  output  file  that  contains  position  and  body  attitude 
formation.  Logical  unit  8  is  an  output  file  that  contains  other  information  useful  to  the  user, 
such  as  the  velocity  angles,  angle  of  attack,  and  acceleration  components.  Each  file  is 
described  in  more  detail  below. 

INPUT  DATA 

Input  for  MOFO  consists  of  two  files,  an  aircraft  characteristics  file  and  a  command 
input  file. 

Aircraft  Characteristics  File 

The  aircraft  characteristics  file  consists  of  both  aerodynamic  data  and  power  data  of 
the  specific  vehicle  being  modeled.  The  aerodynamic  data  consist  of  lift  coefficients,  drag 
coefficients,  and  aerodynamic  time  constants.  The  power  data  consist  of  thrust,  fuel  flow, 
and  thrust-specific  fuel  consumption,  among  other  things.  Aircraft  design  variables  are 
also  required,  such  as  the  wing  surface  area,  base  weight  of  airframe,  and  maximum  fuel 
capacity.  Table  1  shows  how  this  file  is  structured. 

Users  of  the  Advanced  Air-to-Air  System  Performance  Evaluation  Model 
(AASPEM)  may  note  the  similarities  in  the  input  formats  .  Each  line  in  table  1  corresponds 
to  a  record  in  the  file,  except  for  the  representation  of  two-dimensional  tables.  All  variables 
shown  in  the  table  are  explained  in  detail  in  appendix  A  under  the  "Common  Block" 
subsection. 
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Table  1.  Aircraft  input  file  format 


Model  variables 

FORTRAN  format 

blank  line  for  descriptive  label 

NMCHMX,  NAOAMX 

IALT(1)-IALT(NMCHMX) 

SPDM  AX(  1)-SPDM  AX(NMCHMX) 

AERSPD(1  )-AIRSPD(N  AO  AMX) 

AOAMAX(l)-AOAMAX(NAOAMX) 

NALTML,  NMCHML,  SCALE 
blank  line  for  table  desciption 

MCHMILf  1  )-MCHMIL  (NMCHML) 

ALTMIL(1 )  JP  WR  MIL(  1 , 1  )-P  WRMIL(  1  .NMCHML) 

2110 

10X.14I5 

10X.14F5.0 

10X.14F5.0 

10X.14F5.0 

10X,2I5,30X,F10.0 

10X.14F5.0 

ALTMIL(I)  ,PWRMIL(I,1)-PWRMIL(I,NMCHML) 

4X, 15, IX, 1415 

ALTMIL(N  ALTML),PWRMIL(N  ALTML,  1  )-PWRMIL(N  ALTML.NMCHML) 
NALTAB, NMCHAB  .SCALE 
blank  line  for  table  description 

MCHAB(1)-MCHAB(NMCHAB) 

ALT  AB  ( 1 )  ,P  WR  AB  ( 1 , 1  )-P  WR  AB(  1  .NMCHAB) 

10X,2I5,30X,F10.0 

10X.14F5.0 

ALT  AB  (I),PWR  AB(1, 1  )-PWR  AB(I,NMCHAB) 

4X.I5.1X, 1415 

ALTAB(NALTAB),PWRAB(NALTAB,1)-PWRAB(NALTAB,NMCHAB) 
NALTML, NMCHML, SCALE 
blank  line  for  table  descripuon 

MCHMIL(1)-MCHMIL  (NMCHML) 

ALTMIL(l)  ,FULMIL(  1 , 1  )-FULMIL(  1  ,NMC  HML) 

10X,2I5,30X,F10.0 

10X.14F5.0 

ALTMIL(I)>FULMIL(I11)-FULM3L(LNMCHML) 

4X, 15, IX, 1415 

ALTMIL(N  ALTML)  ,FULMIL(N  ALTML ,  1  )-FULMIL(N  ALTML.NMCHML) 
NALTAB, NMCHAB  .SCALE 
blank  line  for  table  description 

MCH  AB  ( 1)-MCHAB  (NMCHAB) 

10X,2I5,30X,F10.0 

10X.14F5.0 
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Table  1.  (Continued) 


Model  variables 

FORTRAN  format 

ALTAB(  1 )  JFULAB  ( 1,1  )-FUL  AB  (1  .NMCHAB) 

ALTAB(I)rFULAB(I,l)-FULAB(I,NMCHAB) 

4X,I5,1X,14I5 

ALTAB(NALTAB)JTXAB(NALTAB,1)-FULAB  (NALTAB  ,NMCHAB) 
NCLCD,NMCHCD, SCALE 
blank  line  for  table  description 

M  ACHCD(  l)-MACHCD(NMCHCD) 
CLFTCD(1),CDRAG(1,1)-CDRAG(1,NMCHCD) 

10X,2I5,30X,F10.0 

10X,14F5.0 

CLFTCD(I),CDRAG(I,1)-CDRAG(I,NMCHCD) 

4X,F5.0,1X,14F5.0 

CLFTCD(NCLCD),CDRAG(NCLCD,1)-CDRAG(NCLCD,NMCHCD) 
NAOACL,NMCHCL, SCALE 
blank  line  for  table  description 

MACHCL(1)-MACHCL(NMCHCL) 

AO  ACL(1 )  ,CLIFT(  1 ,1  )-CLIFT(  1  ,NMCHCL) 

10X,2I5,30X,F10.0 

10X,14F5,0 

AOACL(I),CLIFT(I,l)-CLIFT(I^MCHCL) 

4X,F5.0,1X,14F5.0 

AOACL(NAOACL),CLfFT(NAOACL,l)-CLIFT(NAOACL,NMCHCL) 

NCTRPM 

PCTRPM(1)-PCTRPM(NCTRPM) 

TSFC(1)-TSFC(NCTRPM) 
blank  line  for  variable  description 

DT, DELTA, KAY 

blank  line  for  variable  description 

AO  AT  AU  ,MUT  AU  ,GMT  AU.CFUTAU.VT  AU 
blank  line  for  variable  description 

AREA ,  WT  AC,  WTFLMX 

110 

10X,14F5.0 

10X,14F5.0 

3F10.0 

5F10.0 

3F15.0 
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MOFO  was  developed  specifically  for  modeling  military  turbojet  aircraft.  The  user 
may  find  that  many  of  the  input  parameters  are  easily  accessible  for  fighter  and  attack 
aircraft.  NATOPS  and  TAC  manuals  are  good  initial  starting  points.  These  texts  can  be 
used  to  extract  such  basic  information  as  the  planform  area  (AREA),  base  weight  (WTAC), 
and  maximum  internal  fuel  storage  (WTFLMX).  On  the  other  hand,  aero-dynamic 
coefficients  and  engine-specific  thrust  and  fuel  flow  data  are  harder  to  come  by.  This 
information  is  usually  obtained  from  the  appropriate  manufacturer  or  NAY  AIR.  Finally, 
there  are  variables  that  must  be  adjusted  until  output  performance  matches  the  user's 
expectation  or  validation  data.  Among  this  class  of  input  variables  are  the  integration  time 
step,  command  guidance  time  step,  and  the  proportional  navigation  constant  (see  [1]  for 
details).  Figure  5  provides  the  user  with  sample  aircraft  input  file  for  the  F/A-18  hornet. 

Command  Input  File 

The  command  input  file  consists  of  the  aircraft's  initial  position  followed  by  a  series 
of  instructions.  These  instructions  constitute  the  commanded  route  on  which  the  trajectory 
is  generated.  Such  commands  as  what  position  to  fly  to  or  what  heading  to  follow  and 
with  what  climb  or  dive  angle  provide  the  user  flexibility  and  control  in  prescribing  a 
desired  behavior.  Users  familiar  with  BLUEMAX II  may  notice  similarities  between  this 
file  and  its  flight  control  file  used  in  the  automated  mode.  The  format  for  this  input  file  is 
provided  in  table  2. 

The  file  is  composed  of  three  parts.  The  first  part  describes  the  initial  position,  the 
direction  and  magnitude  of  the  velocity  vector,  the  angle  of  attack,  and  weight  considera¬ 
tions.  The  first  line  provides  the  user  with  a  line  for  a  flight  path  description.  A  blank  line 
then  follows.  Next,  the  initial  Cartesian  position  of  the  aircraft  and  its  speed  are  provided. 
The  next  record  has  the  velocity  angles  and  the  angle  of  attack.  The  final  record  contains 
the  initial  time  of  the  simulation  and  the  fuel  and  armament  weights. 

The  second  part  consists  of  a  series  of  commands  that  maneuver  the  aircraft. 
Currently,  MOFO  has  one  instruction  named  WAYPT.  WAYPT  flies  the  aircraft  to  the 
next  user-specified  position  along  a  path  that  is  aerodynamically  feasible.  The  user  also 
specifies  the  speed  to  be  attained  when  the  waypoint  is  achieved  and  a  g  limit  to  impose 
constraints  on  directional  changes.  Future  commands  to  be  implemented  are  discussed  in 
detail  in  the  final  section. 

The  third  and  final  part  of  this  input  file  is  the  END  command.  It  concludes  the 
flight  path  generation  and  stops  program  execution.  Figure  6  shows  a  sample  command 
input  file. 

OUTPUT  DATA 

Output  from  the  model  is  contained  in  two  files.  One  file  (named  POSITION  .OUT 
in  the  sample  ".COM"  file)  provides  the  user  with  simulation  time,  aircraft  positions  (X,  Y, 
and  Z),  speed,  and  body  attitude  (*b,  V  and  M-b)-  Figure  7  is  a  sample  of  this  output  using 
the  format  (IX,  5F11.2,  3F11.4). 

The  second  output  file  (named  VELOCITY. OUT  in  the  sample  ".COM”  file) 
provides  the  user  with  other  useful  information  on  the  aircraft.  This  file  contains  weight, 
velocity  angles  (%,  y,  and  p.),  angle  of  attack  (a),  thrust,  and  the  acceleration  components 
(longitudinal,  pitch,  yaw,  and  total  accelerations).  Figure  8  is  a  sample  of  this  output  using 
the  format  (IX,  1 1F10.2). 
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F/A-18  AIRCRAFT  CHARACTERISTIC  FILE 


12 

12 

0 

50001 

0000 

150002000025000 

300003500040000450005000055000 

1 . 1 

1.24 

1.38 

1.51 

1.65 

1.74 

1.87 

1  .  98 

1 . 98 

1.99  2.0  2.0 

0.2 

0.6 

0.8 

0.85 

0.9 

0.95 

1 .05 

1.1 

1.2 

1.4  1.6  1.8 

36.4 

35.5 

35.8 

35.4 

35 . 0 

32.0 

21 . 2 

19.8 

16 . 8 

14.6  14.8  14.6 

12 

10 

2.0 

THRUST, 

MILITARY  : 

POWER 

POUNDS 

0. 

0.2 

0  .  4 

0.6 

0.8 

1 . 0 

1.2 

1.4 

1.6 

1.8 

0 

9556 

9541 

9693 

9549 

9370 

7895 

4259 

523- 

1576- 

6745 

5000 

8255 

8233 

8646 

8723 

8868 

8585 

5772 

2761 

496- 

■3620 

10000 

7304 

7322 

7704 

7788 

8031 

7974 

6947 

4439 

2362 

-77 

15000 

6111 

6140 

6497 

6845 

7159 

7313 

6648 

5411 

3089 

1001 

20000 

5025 

5057 

5410 

5729 

6265 

6514 

6287 

5888 

4212 

1719 

25000 

4035 

4085 

4415 

4738 

5297 

5694 

5740 

5680 

4781 

2262 

30000 

3213 

3230 

3537 

3849 

4365 

4907 

5081 

5372 

5124 

3327 

35000 

2559 

2575 

2795 

3050 

3548 

4014 

4363 

4816 

4834 

3808 

40000 

2010 

2001 

2145 

2351 

2746 

3131 

3457 

3847 

3906 

3195 

45000 

1625 

1619 

1663 

1781 

2096 

2403 

2679 

3000 

3045 

2486 

50000 

1284 

1288 

1331 

1404 

1585 

1826 

2070 

2306 

2358 

1912 

55000 

951 

970 

1031 

1115 

1253 

1399 

1562 

1770 

1816 

1464 

12 

10 

2.0 

THRUST , 

AFTER  BURNER 

POUNDS 

0. 

0.2 

0 . 4 

0.6 

0.8 

1.0 

1.2 

1.4 

1.6 

1.8 

0  143981489815721163101716516431145031253710748-4118 


5000 

12423128391398314760159281697015752148371322310087 

10000 

10943113251234813095142451541516584155611419013233 

15000 

9171 

61401041511417125801386115132155331361412762 

20000 

5750 

208 

5854  9544109231219713678149581379411996 

25000 

6143 

6410 

79  -99 

91661057712101135031327511174 

30000 

4954 

5139 

5818  4204 

119 

896110472120621262711278 

35000 

3991 

4147 

4630  5207 

3856 

7306 

88301 

04181118210695 

40000 

3138 

3227 

413  4022 

385 

3429 

6978 

8278 

8945 

8692 

45000 

2487 

2585 

2757  397 

2286 

4425 

5425 

6441 

6966 

6762 

50000 

1973 

2047 

2188  2391 

2789 

3376 

4182 

4963 

5394 

5245 

55000 

1438 

1524 

1677  1876 

2176 

2569 

3165 

3803 

4086 

4047 

12 

10 

2.0 

FUEL  FLOW,  MILITARY  POWER 

POUNDS/HOUR 

0. 

0.2 

0.4  0.6 

0.8 

1.0 

1.2 

1.4 

1.6 

1.8 

0 

7963 

8523 

9310  9929 

1049210277 

9945 

9480 

9123 

8740 

5000 

6861 

7345 

8189  8894 

97311039710201 

9705 

9332 

8978 

10000 

6081 

6509 

7243  7804 

8595 

9345 

10105 

9696 

9015 

9174 

15000 

5053 

5414 

6042  6779 

7495 

8313 

9069 

9347 

8264 

8300 

20000 

4145 

4440 

4995  5611 

6462 

7239 

8085 

8708 

8082 

7420 

25000 

3299 

3559 

4045  4597 

5386 

6218 

7087 

7750 

7668 

6633 

30000 

2593 

2786 

3196  3699 

4380 

5261 

6099 

6841 

7107 

6486 

35000 

2029 

2182 

2495  2891 

3527 

4234 

5155 

5902 

6255 

6105 

40000 

1599 

1697 

1924  2238 

2739 

3323 

4098 

4717 

5029 

4970 

45000 

1318 

1396 

1521  1722 

2118 

2583 

3207 

3698 

3947 

3900 

50000 

1063 

1132 

1241  1392 

1635 

2000 

2498 

2885 

3094 

3056 

55000 

801 

869 

980  1124 

1326 

1563 

1925 

2251 

2421 

2391 

12 

10 

2.0 

FUEL  FLOW,  AFTER  BURNER 

POUNDS/HOUR 

0  . 

0.2 

0.4  0.6 

0.8 

1 . 0 

1.2 

1.4 

1.6 

1.8 

0  285753069833770367984042641966440344402943694  8881 
5000  24251260082922732243363524059042943446374434644116 
10000  21008225162527127777314123552940707426154325044678 
15000  17445  54142094023673268683085035350391833880140496 
20000  8298  607  645919442227272626830637351263585335725 

25000  1175312601  587  574187102210826173301753228631470 


Figure  5.  Sample  aircraft  characteristics  input  file 
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30000  95801021511671  4925  5881835822029258722870728944 

35000  7806  8322  937410577  46821467718230217242430825735 

40000  6248  6584  635  8273  623  442714463173001945420772 

45000  5075  5395  5809  643  2969  905611353136031530716334 

50000  4128  4379  4726  5168  5894  7065  8889106391203812831 

55000  3184  3413  3757  4178  4749  5549  6888  8340  899110072 

12  11  1.0 

DRAG  COEFFICIENTS 

0.0  0.85  0.9  0.95  1.0  1.05  1.2  1.4  1.6  1.8  2.0 

0.0  .0276.0237.0264.0368.0510.0531.0539.0563.0571.0580.0587 
0.05  .0272.0234.0261.0367.0500.0519.0525.0543.0548.0551.0552 
0.1  .0276.0237.0264.0368.0502.0519.0526.0546.0552.0557.0558 
0.2  .0300.0261.0287.0388.0528.0552.0567.0603.0628.0657.0684 
0.3  .0348.0309.0334.0445.0613.0625.0667.0733.0796.0861.0928 
0.4  .0424.0386.0416.0557.0742.0751.0812.0908.1016.1136.1240 
0.5  .0526.0499.0534.0712.0922.0941.1012.1158.1321.1486.1666 
0.6  .0662.0651.0704.0927.1147.1171.1277.1493.1716.1946.2175 
0.8  .1028.1101.1274.1522.1742.1771.1992.2393.2801.2871.3245 
1.0  .1466.1984.2184.2407.2602.2646.3032.3763.2801.2871.3245 
1.2  .1202.3264.3454.3697.3932.4021.4497.4143.2801.2871.3245 
1.35  .3502.4524.4724.4987.5242.5356.4947.4143.2801.2871.3245 
11  12  1.0 
LIFT  COEFFICIENTS 


0.2 

0.6 

0.8 

0.85 

0.9 

0.95 

1.05 

1.1 

1.2 

1.4 

1.6 

1.8 

0.0 

-.01 

-.02 

-.03 

-.02 

-.05 

-.04 

0.00 

0.00 

0.00 

0.00 

.01 

.01 

4.0 

.33 

.  39 

.44 

.47 

.51 

.51 

.43 

.34 

.  32 

.28 

.24 

.22 

8.0 

.69 

.81 

.85 

.89 

.92 

.88 

.74 

.61 

.  61 

.  53 

.  44 

.41 

12.0 

1.04 

1 .15 

1 .  08 

1 .05 

1.08 

1.18 

1.00 

.86 

.87 

.75 

.63 

.60 

16.0 

1.22 

1.29 

1.18 

1.17 

1.19 

1.35 

1.21 

1.10 

1.10 

.95 

.80 

.77 

20.0 

1.36 

1 .  42 

1.29 

1.28 

1 . 30 

1.46 

1.35 

1.31 

1.30 

.95 

.80 

.77 

24.0 

1.49 

1.54 

1.43 

1  .  39 

1 .40 

1.53 

1.48 

1.31 

1.30 

.95 

.80 

.77 

28 . 0 

1.61 

1.66 

1.55 

1.52 

1.53 

1 .  57 

1.48 

1.31 

1 .  30 

.95 

.80 

.77 

32.0 

1.76 

1.78 

1.63 

1.61 

1.60 

1.58 

1 .46 

1.31 

1.30 

.95 

.80 

.77 

36  .  0 

1.82 

1.81 

1.69 

1.64 

1.61 

1.57 

1 . 48 

1 .  31 

1.30 

.95 

.80 

.77 

40.0 

8 

1.80 

1.76 

1.62 

1.60 

1.59 

1.57 

1.48 

1.31 

1.30 

.95 

.80 

.77 

.  325 

.  4 

.  5 

.6 

.7 

.8 

.9 

1.0 

2.75 

2 . 4 

1.95 

1.6 

1.3 

1 . 1 

.  99 

1 . 0 

DT(sec)  DELTA(sec)  KAY 
0.1000  1.0000  16.0000 

AOATAU ( sec )  MUTAU (sec)  GMTAU(sec)  CHITAU(sec)  VTAU(sec) 
.10000  .1000  1.0000  1.0000  1.0000 

PLANFORM  AREA(sqft)  WTAC(lbs)  WTFLMX(lbs) 

400.0000  23500.00  10810.0 

Figure  5.  (Continued) 
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Table  2.  Command  file  format 


Model  variables 

Units 

FORTRAN  format 

blank  line  for  profile  label 

blank  line 

— 

— 

X0,Y0,Z0,V0 

m,m,m,  m/sec 

4F10.0 

chio,gmo,muo,alphao 

deg, deg, deg, deg 

4F10.0 

TO, WTFUEL, BINGO, WTARM 

sec, lb, lb, lb 

4F10.0 

blank  line 

— 

— 

WAYPT 

command  name  only 

A 

XF,YF,ZF,VF,GEES 

m,m,m, m/sec, — 

5F10.0 

blank  line 

— 

— 

END 

command  to  end 

A 

SAMPLE  RUN  FOR  F/A-18 


-10000.0 

90.0 

0 . 0 

0  .  0 
0.0 

20000 . 0 

172 .24 
0.0 
0.0 

139 . 01 

4.0 

0.0 

WAYPOINT 

0.0 

0 . 0 

172.24 

139.01 

5.0 

WAYPOINT 

0  .  0 

10000.0 

172.24 

139 . 01 

2 . 0 

WAYPOINT 

10000.0 

10000.0 

172 .24 

139.01 

5.0 

END 

Figure  6.  Sample  command  input  file 
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X 

(  m  )  , 

.  Y  ( in  )  ,  Z  (  m  )  , 

-  V  ( io/s  )  , 

BODY  YAW 

( deg . ) ,  BODY 

PITCH  (deg. 

) ,  BODY 

RO  L  L  (deg. 

0  . 

.  0  0 

-10000  .  00 

0.00 

172.24 

139.01 

90.0000 

4.0000 

0.0000 

1 . 

.  00 

-9860.99 

0.00 

172.26 

139.01 

90.0000 

3.8948 

0  .  00  00 

2  . 

.  0  0 

-972  1  .98 

0.00 

172.27 

139.01 

90.0000 

3.8927 

0.0000 

3  . 

.  0  0 

-9582.97 

0.00 

172.29 

13  9  .01 

90.0000 

3.8911 

0.0000 

4  . 

.  0  0 

-9443.96 

0.00 

172.30 

139  .01 

90.0000 

3  .  8  8  9  7 

0.0000 

5  . 

.  0  0 

-9304.95 

0.00 

17  2.30 

139.01 

90.0000 

3  .  8  8  8  7 

0.0000 

6  . 

.  00 

-9165.94 

0.00 

172.3  1 

139.01 

90.0000 

3  .  8  879 

0.0000 

7 

.  0  0 

-9026.93 

0.00 

172.3  1 

139.01 

90.0000 

3.8872 

0  .0  000 

8.00 

-8887.92 

0.00 

172.3  2 

139.01 

90.0000 

3  .  8  8  6  7 

0.0000 

9 

.  00 

-8748.91 

0.00 

172.32 

139.01 

90.0000 

3.8864 

0.0000 

1  0 

.  00 

-8609.90 

0.00 

172.32 

139.01 

90.0000 

3.8860 

0.0000 

1 1  . 

.  00 

-8470.89 

0.00 

17  2.32 

139.01 

90.0000 

3  .  8  8  5  8 

0.0000 

12.00 

-8  3  3  1  .88 

0.00 

172.32 

139.01 

90.0000 

3.8  8  5  6 

0.0000 

13.00 

-8192.87 

0.00 

172.32 

139  .01 

90.0000 

3.8854 

0.0000 

1  4  . 

.  00 

-8053.86 

0.00 

17  2.32 

139  .01 

90.0000 

3.8  8  5  2 

0.0000 

15.00 

-7914.85 

0.00 

172.31 

139.01 

90.0000 

3.8850 

0.0000 

16  . 

.  0  0 

-7775 . 84 

0.00 

172.3  1 

139.01 

90.0000 

3.8848 

0.0000 

17 

.  0  0 

-7636.83 

0.00 

172.31 

139.01 

9  0  .0  0  0  0 

3.8847 

0.0000 

1  8 

.  0  0 

-74  9  7  .  82 

0.00 

172.31 

139.01 

90.0000 

3.8845 

0.0000 

1  9 

.  00 

-7358.81 

0.00 

172.31 

139.01 

90.0000 

3  .  8  8  4  4 

0.0000 

20.00 

-7219.80 

0.00 

17  2.31 

139.01 

90.0000 

3.8842 

0.0000 

21 

.  0  0 

-70  8  0  .  79 

0.00 

172.31 

139.01 

90.0000 

3  .  8  8  4  1 

0.0000 

22 

.  00 

-6941.78 

0.00 

172  .  31 

139.01 

90.0000 

3  .  8  8  3  9 

0.0000 

23 

.  00 

-6802.77 

0.00 

172.31 

139.01 

90.0000 

3  .  8  8  3  7 

0.0000 

24 

.  00 

-6663.76 
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Figure  7.  Sample  position  and  attitude  output 
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Figure  7.  (Continued) 
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Figure  7.  (Continued) 
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Figure  8.  Sample  velocity  output 
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Figure  8.  (Continued) 


CONCLUSIONS 


This  section  briefly  describes  improvements  to  MOFO,  certain  tradeoffs  between 
fidelity  and  run  time,  and  future  work.  The  following  information  supplements  that  in  [3]. 

VARIABLE  RANGE 

At  present,  MOFO  computes  the  trajectory  of  an  aircraft  while  en  route  toward  a 
specific  waypoint  in  space.  It  determines  that  the  waypoint  has  been  achieved  if  the 
trajectory  comes  within  a  predetermined  "basket."  This  basket  is  described  as  a  set  range 
from  the  waypoint.  If  the  user  either  wishes  to  have  more  flexibility  in  this  predetermined 
range  or  to  have  the  aircraft  come  as  close  as  possible  to  the  waypoint,  two  suggestions  are 
given. 


If  more  flexibility  in  controlling  the  course  of  the  trajectory  is  desired,  a  simple 
modification  should  suffice.  For  example,  the  user  may  input  a  specific  range  at  the 
beginning  of  the  command  file  or  provide  a  range  to  accompany  each  WAYPT  command. 

Under  some  circumstances,  however,  it  may  be  desirable  to  have  the  aircraft  follow 
the  desired  trajectory  more  closely.  In  this  instance,  the  basket  essentially  shrinks  to  a 
point.  This  situation  is  more  difficult  to  model.  This  capability  would  depend  on  the 
severity  of  the  maneuver  and  the  current  values  of  the  velocity,  DELTA,  and  DT.  The 
mathematical  development  is  left  for  future  users. 

OTHER  CONTROL  OPTIONS 

Currently,  MOFO  has  only  one  control  option,  WAYPT.  Users  may  wish  to  have 
more  flexibility  in  controlling  the  attitude  and  direction  of  the  aircraft.  The  source  code  of 
subroutine  CONTROL  in  appendix  B  contains  other  control  routines  that  have  been 
"commented"  out.  These  routines  are  yet  to  be  developed  but  provide  a  foundation  for 
needed  improvements.  They  are: 

•  HEADIN.  This  subroutine  would  allow  control  of  the  velocity  heading 
angle,  %. 

•  PITCH.  This  subroutine  would  allow  control  of  the  velocity  pitch  angle,  y. 

•  HDGPTC.  This  subroutine  would  allow  control  of  both  the  velocity 
heading  and  velocity  pitch  angle,  %  and  y,  respectively. 

•  DRPSTR.  This  subroutine  would  allow  control  of  the  dropping  of  external 
stores  from  the  aircraft  amounting  to  a  weight  reduction. 

•  TERFOL.  This  subroutine  would  allow  the  aircraft  to  perform  terrain 
following  to  a  waypoint.  It  would  require  the  acquisition  and  manipulation 
of  terrain  elevation  data. 

Another  control  routine  not  included  in  the  source  code  will  also  need  to  be 
developed  if  necessity  dictates.  Namely,  MOFO  reads  in  thrust  and  fuel  flow  data  for  both 
military  and  afterburner  power.  However,  the  model  computes  thrust  under  military  power 
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only.  A  switch  would  need  to  be  implemented  if  the  user  wishes  to  fly  the  aircraft  in 
afterburner  power. 

FIDELITY  VERSUS  CPU  TIME 

As  with  all  models,  MOFO  is  not  without  the  tradeoff  that  exists  between  accuracy 
and  run  time.  In  MOFO,  a  reasonable  compromise  between  fidelity  and  simplicity  has  been 
achieved  in  bringing  the  model  to  its  current  one-to-one  ratio  of  actual-to-simulated  time. 
For  instance,  MOFO  does  not  model  the  aerodynamic  response  at  the  control  surface  level 
(simulating  the  response  to  deflections  of  the  rudder,  fins,  ailerons,  etc.,  as  described  in 

[4] ).  This  level  of  detail  was  not  required  for  the  intended  use  of  the  model. 

If  a  shorter  run  time  is  required  with  position  and  attitude  information  at  a  lesser 
degree  of  accuracy,  some  suggestions  are  in  order.  One  approach  would  be  to  convert  the 
double  precision  format  into  single  precision.  Another  would  be  to  use  linear  interpolation 
as  opposed  to  cubic  and  bicubic  spline  interpolation.  A  third  possibility  would  be  an 
integration  algorithm  other  than  Runge-Kutta.  Although  the  Runge-Kutta  method  is  self¬ 
starting,  a  great  deal  of  computer  time  is  required  for  the  four  intermediate  function 
evaluation  steps.  An  alternative  may  be  to  use  the  Adams  formulas  (such  as  those  found  in 

[5] )  after  the  first  few  Runge-Kutta  iterations. 

CURRENT  WORK 

At  the  time  of  this  writing,  work  is  under  way  both  to  validate  and  modify  MOFO. 
The  validation  effort  has  involved  the  comparison  of  performance  predictions  from  MOFO 
and  the  generally  accepted  flight  path  generator,  BLUEMAX  II,  to  aerodynamic  charts 
found  in  N  ATOPS  manuals.  MOFO  is  being  modified  to  include  implementing  a  variable 
DELTA  time  step  that  depends  on  the  maneuvers  being  performed  by  the  aircraft  (see  [6]). 
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MODULE  AND  VARIABLE  DESCRIPTIONS 


PROGRAM  MOFO 

The  program  MOFO  is  the  main  module  of  the  MOFO  model.  MOFO  begins  by 
initializing  several  variables  that  are  used  throughout  the  rest  of  the  program.  It  then  reads 
in  the  aircraft  input  data  file.  This  file  consists  of  one-dimensional  and  two-dimensional 
tables.  Program  MOFO  reads  in  the  one-dimensional  tables  and  calls  subroutine  TABIN  to 
read  the  two-dimensional  tables.  Subroutines  SPLINE  and  SPLIE2  are  called  concurrently 
to  create  second  derivative  tables  that  will  be  used  later  when  cubic  and  bicubic  spline 
interpolation  is  required.  For  some  of  the  tables,  subroutine  PARDER  is  called  to  create 
partial  derivative  tables  to  be  used  later  in  the  root-finding  subroutine,  NWTN1D.  Finally, 
program  MOFO  calls  subroutine  CONTROL  to  begin  processing  the  command  file. 

MOFO  Common  Variables 

The  following  are  variables  placed  in  common  blocks  that  are  used  throughout  the 
model  MOFO.  These  common  blocks  are  located  in  the  file  called  CMNBLK.CMN.  Most 
routines  have  an  "INCLUDE"  statement  to  incorporate  this  structure. 

Parameters 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 

Common  Block  I AERO/ 

AIRSPD(NMAX) — array  of  Mach  numbers  of  length  NAOAMX 

AOAMAX(NMAX) — array  of  the  maximum  angle  of  attack  for  each  Mach 
number  in  AIRSPD(NMAX)  (degrees) 

AAMXSD(NMAX) — array  that  contains  the  second  derivatives  of 
AOAMAX(NMAX)  at  each  Mach  number  in  AIRSPD(NMAX) 

CLFTCD(NMAX) — array  of  lift  coefficients  of  length  NCLCD 

MACHCD(NMAX) — array  of  Mach  numbers  of  length  NMCHCD 

CD — interpolated  value  of  the  drag  coefficient,  Cd 

CDRAG(NMAX,NMAX) — two-dimensional  array  of  the  drag  coefficients 
as  a  function  of  CLFTCD(NMAX)  and  MACHCD(NMAX),  of  size 
NCLCD  by  NMCHCD 

CDRGSD(NMAX,NMAX) — two-dimensional  array  that  contains  the  second 
derivatives  of  the  drag  coefficients  associated  with  each  value  in 
CDGRAG(NMAX,NMAX) 

AOACL(NMAX) — array  of  angles  of  attack  of  length  NAOACL  (degrees) 

MACHCL(NMAX) — array  of  Mach  numbers  of  length  NMCHCL 
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CL — interpolated  value  of  the  lift  coefficient,  Cl 

CLIFT(NMAX,NMAX) — two-dimensional  array  of  the  lift  coefficients  as  a 
function  of  AOACL(NMAX)  and  MACHCL(NMAX),  of  size  N  AO  ACL 
by  NMCHCL 

CLFTSD(NMAX,NMAX) — two-dimensional  array  that  contains  the  second 
derivatives  of  the  lift  coefficients  associated  with  each  value  of 
CLIFT(NMAX,NMAX) 

DERCD — interpolated  value  of  the  partial  derivative,  0Cd/0Cl 

DCDDCL(NMAX,NMAX) — two-dimensional  array  of  the  partial  derivative, 
0Qd/3Cl,  of  size  NCLCD  and  NMCHCD  associated  with  each  value  of 
CDRAG(NMAX,NMAX) 

DCDDCLSD(NMAX,NMAX) — two-dimensional  array  that  contains  the 
second  derivatives  of  each  value  in  DCDDCL(NMAX,NMAX) 

DERCL — interpolated  value  of  the  partial  derivative,  dCjJda 

DCLDA(NMAX,NMAX) — two-dimensional  array  of  the  partial  derivative, 
3Cl/9cx,  associated  with  each  value  of  CL1FT(NMAX,NMAX) 

DCLDASD(NMAX,NMAX)-two-dimensional  array  that  contains  the  second 
derivatives  of  each  value  in  DCLDA(NMAX,NMAX) 

Common  Block  IAIRCFTI 

ACM — total  mass  of  aircraft,  fuel,  and  armaments  (slugs) 

AREA — planform  area  (ft2) 

BINGO — amount  of  fuel  left  that  signals  the  pilot  to  return  to  base  (lb) 

WTAC — base  weight  of  the  aircraft  (lb) 

WTFLMX — maximum  weight  of  fuel  that  can  be  loaded  on  the  planform  (lb) 

WTFUEL — initial  weight  of  fuel  placed  in  the  command  file  (lb) 

Common  Block  I  AO  Al 

ALPHA — present  angle  of  attack  after  time  filtering 

ALPHA 0 — initial  or  previous  angle  of  attack 

AOATAU — aerodynamic  time  constant  for  angle  of  attack  filtering  (seconds) 
Common  Block  /CONSTS/ 

G — acceleration  due  to  gravity  (32.17  ft/sec2) 
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COUNT — value  of  1/DT 
ICOUNT — integer  value  of  COUNT 

MTOFT — a  conversion  factor:  meters  to  feet  (3.280840  ft/m) 

P0 — reference  pressure  at  sea  level  (21 16.2  lb/ft2) 

PI— k  (3.14159265358979) 

RHOO — density  of  atmosphere  at  sea  level  (2.377  x  10-3  slugs/fO) 

ROOT — value  of  1/3.5  used  in  the  RPM  calculation  in  routine  AERO 

RTODEG — a  conversion  factor:  radians  to  degrees  (180/7t  degrees/radian) 
Common  Block  ICONTRL/ 

CHIDOT — present  commanded  heading  rate  after  time  filtering 

CHIDOTO — previous  commanded  heading  rate 

CHITAU — time  constant  for  filtering  yaw  rate  commands  (seconds) 

GMDOT — present  commanded  pitch  rate  after  time  filtering 
GMDOTO — previous  commanded  pitch  rate 

GMTAU — time  constant  for  filtering  pitch  rate  commands  (seconds) 

VDOT — present  commanded  acceleration  after  time  filtering 
VDOTO — previous  commanded  acceleration 

VTAU — time  constant  for  filtering  longitudinal  accleration  commands 
(seconds) 

DELTA — control  time  step  interval  between  acceleration  commands  (seconds) 
DT — differential  time  step  for  the  Runge-Kutta  algorithm  (seconds) 

GEES — magnitude  of  the  load  factor  that  the  platform  experiences 
KAY — constant  for  proportional  navigation 
Common  Block  /CYCLE! 

MV — momentum  (mass  times  acceleration  due  to  velocity) 

MG — weight  (mass  times  gravity) 
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REFOR — reference  force  (,5*p*V2*AREA) 
SNG — sin(y);  sine  of  the  velocity  pitch  angle 
CSG — cosCy);  cosine  of  the  velocity  pitch  angle 


SNCHI — sin  (%);  sine  of  the  velocity  heading  angle 

CSCHI — cos(x);  cosine  of  the  velocity  heading  angle 
CT — interpolated  value  of  the  thrust  specific  fuel  consumption 
Common  Block  IFNCTNI 

FI — X;  the  velocity  component  in  the  X  direction 

F2 — Y;  the  velocity  component  in  the  Y  direction 

F3 — Z;  the  velocity  component  in  the  Z  direction 

F A — V;  the  longitudinal  acceleration  component  along  the  velocity  vector 

F5 — J,  the  pitch  rate 

F6 — X  the  yaw  rate 

F7 — m;  the  rate  of  change  of  the  aircraft's  mass 
Common  Block  /POWER/ 

ALT(NMAX) — array  of  altitudes  of  length  NMCHMX  (feet) 

SPDMAX(NMAX) — array  of  maximum  Mach  numbers  for  each  altitude  in 
ALT(NMAX) 

SPMXSD(NMAX) — array  of  second  derivatives  associated  with  each  value  in 
SPDMAX(NMAX) 

ALTMIL(NMAX) — array  of  altitudes  of  length  NALTML  (feet) 
MCHMIL(NMAX) — array  of  Mach  numbers  of  length  NMCHML 
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PWRMIL(NMAX,NMAX) — two-dimensional  array  of  maximum  thrust 
under  military  power  as  a  function  of  ALTMIL(NMAX)  and 
MCHMIL(NMAX)  of  size  NALTML  by  NMCHML  (lb) 

PRMLSD(NMAX,NMAX) — two-dimensional  array  of  the  second  derivatives 
associated  with  each  value  in  PWRMIL(NMAX,NMAX) 

FULMIL(NMAX,NMAX) — two-dimensional  array  of  maximum  fuel  flow 
rates  under  military  power  as  a  function  of  ALTMEL(NMAX)  and 
MCHMIL(NMAX)  of  size  NALTML  by  NMCHML  (lb/hr) 

FLMLSD(NMAX,NMAX) — two-dimensional  array  of  the  second  derivatives 
associated  with  each  value  in  FULMIL(NMAX,NMAX) 

ALTAB(NMAX) — array  of  altitudes  of  length  NALTAB  (feet) 

MCHAB(NMAX) — array  of  Mach  numbers  of  length  NMCHAB 

PWRAB(NMAX,NMAX) — two-dimensional  array  of  maximum  thrust  under 
afterburner  power  as  a  function  of  ALTAB(NMAX)  and 
MCHAB(NMAX)  of  size  NALTAB  by  NMCHAB  (lb) 

PRABSD(NMAX,NMAX) — two-dimensional  array  of  the  second  derivatives 
associated  with  each  value  in  PWRAB(NMAX,NMAX) 

FULAB(NMAX,NMAX) — two-dimensional  array  of  fuel  flow  rates  under 
afterburner  power  as  a  function  of  ALTAB(NMAX)  and 
MCHAB(NMAX)  of  size  NALTAB  by  NMCHAB  (Ib/hr) 

FLABSD(NMAX,NMAX) — two-dimensional  array  of  the  second  derivatives 
associated  with  each  value  in  FULAB(NMAX,NMAX) 

PCTRPM(NMAX) — array  of  percentages  of  maximum  revolutions  per  minute 
(RPMs)  of  the  turbojet  engine  of  length  NCTRPM 

TSFC(NMAX) — array  of  thrust  specific  fuel  consumption  for  each  value  in 
PCTRPM(NMAX) 

TSFCSD(NMAX) — array  of  the  second  derivatives  associated  with  each 
value  in  TSFC(NMAX)  (lb/hr/lb) 

MACH — Mach  number  at  the  aircraft's  present  altitude 

TH — computed  thrust  value  needed  to  satisfy  the  equations  of  motion 

THMAX — interpolated  value  of  the  maximum  thrust  at  the  aircraft's  present 
altitude  and  Mach  number 
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Common  Block  /ROLL/ 


MU — roll  angle  after  time  filtering 
MUO — initial  or  previous  roll  angle 

MUTAU — aerodynamic  time  constant  for  roll  angle  filtering  (seconds) 
Common  Block  ITBLMAXI 

NMCHMX— the  number  of  values  in  ALT(NMAX),  SPDMAX(NMAX), 
and  SPMXSD(NMAX) 

NAOAMX — the  number  of  values  in  AIRSPD(NMAX),  AOAMAX(NMAX), 
and  AAMXSD(NMAX) 

NALTML— the  number  of  values  in  ALTMEL(NMAX) 

NMCHML— the  number  of  values  in  MCHMIL(NMAX) 

NALTAB  —the  number  of  values  in  ALTAB(NMAX) 

NMCHAB — the  number  of  values  in  MCHAB(NMAX) 

NCLCD — the  number  of  values  in  CLFTCD(NMAX) 

NMCHCD — the  number  of  values  in  MACHCD(NMAX) 

NAOACL — the  number  of  values  in  AOACL(NMAX) 

NMCHCL — the  number  of  values  in  MACHCL(NMAX) 

NCTRPM— the  number  of  values  in  PCTRPM(NMAX),  TSFC(NMAX),  and 
TSFCSD(NMAX) 

MOFO  Local  Variables 

The  following  variables  are  intrinsic  only  to  program  MOFO: 

I — counter  for  DO  loop 

IALT(NMAX) — array  of  altitudes  in  integer  form  that  is  later  converted  to 
double  precision  real  and  placed  in  array  ALT(NMAX)  (feet). 

SUBROUTINE  AERO 

The  subroutine  AERO  is  called  by  the  various  control  routines.  Its  function  is  to 
compute  the  aerodynamic  forces  based  on  the  commanded  accelerations.  Required  values 
of  the  angle  of  attack,  the  velocity  roll  angle,  and  thrust  are  determined  via  a  root-finding 
algorithm  in  subroutine  NWTN1D.  The  passed  parameters  are: 

X — the  present  value  of  the  X  coordinate 

Y — the  present  value  of  the  Y  coordinate 
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Z — the  present  value  of  the  Z  coordinate 
V — the  present  magnitude  of  the  velocity  vector 

GM — the  present  velocity  pitch  angle,  y 

CHI — the  present  velocity  yaw  angle,  % 

T — the  present  time  on  the  model’s  clock. 

AERO's  local  variables  are: 

H — the  altitude  of  the  aircraft 

ROH — the  atmospheric  density,  p,  at  altitude  H 

VSQ — the  square  of  the  present  speed,  V2 

MUI — the  instructed  roll  angle,  pi,  calculated  from  the  directional 
accelerations,  CHIDOT  and  GMDOT 

ALPHAI — the  instructed  angle  of  attack,  oq,  returned  from  subroutine 
NWTN1D;  the  present  ALPHA  is  used  as  an  initial  guess 

TEMALF — a  temporary  place  holder  when  limiting  ALPHAI  between  -%  and 

K 

ALFMAX — the  maximum  angle  of  attack  at  the  present  Mach  number 

TH — the  present  thrust  needed  to  satisfy  the  commanded  aerodynamic 
conditions 

THMAX — the  maximum  thrust  at  the  present  altitude  and  Mach  number 

RPM— the  percentage  of  maximum  revolutions  per  minute  (RPMs)  of  the 
turbojet  engine  as  a  function  of  TH  and  THMAX. 

SUBROUTINE  ATTITUDE 

The  subroutine  ATTITUDE  is  called  by  subroutine  CONTROL  and  subroutine  RK. 
At  present,  this  subroutine  is  called  once  every  second  of  simulation  time.  It  determines  the 
body  attitude  of  the  aircraft  and  prints  the  results  into  a  file.  The  passed  parameters  are: 

XI — simulation  time 

X2 — X 

X3 — Y 
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X4 — Z 


X5 — V 
X6 — y 
X7 — % 

X8 — a 
X9 — ji. 

Attitude's  local  variables  are: 

CAL — cosine  of  a 
SAL — sine  of  a 
CGM — cosine  of  y 
SGM — sine  of  y 
CCH — cosine  of  % 

SCH — sine  of  % 

CMU — cosine  of  q 
SMU — sine  of  q 

AA — X-component  of  the  %  unit  vector  in  the  internal  local  horizon 
system 

BB — Y-component  of  the  ib  unit  vector  in  the  internal  local  horizon  system 
CC — Z-component  of  the  ib  unit  vector  in  the  internal  local  horizon  system 
BCHI — the  body  yaw  angle  with  respect  to  the  local  horizon 
BGM — the  body  pitch  angle  with  respect  to  the  local  horizon. 

SUBROUTINE  CONTROL 

The  subroutine  CONTROL  is  called  from  program  MOFO.  Its  task  it  to  process  the 
command  file.  First,  the  initial  conditions  of  the  aircraft  are  read  in.  Then  subsequent 


A-8 


orders  are  read  in  and  transferred  to  the  appropriate  subroutine  to  complete  the  command. 
The  subroutine's  local  variables  are: 

XO — initial  X  position  of  aircraft  (m) 

YO — initial  Y  position  of  aircraft  (m) 

ZO — initial  Z  position  of  aircraft  (m) 

VO — initial  speed  of  aircraft  (m/sec) 

GMO — initial  velocity  pitch  angle  of  aircraft  (degrees) 

CHIO — initial  velocity  yaw  angle  of  aircraft  (degrees) 

TO — starting  time  (seconds) 

WTARM — weight  of  external  stores  carried  on  aircraft  (lb) 

X — XO  value  in  the  internal  coordinate  system  (feet) 

Y — YO  value  in  the  internal  coordinate  system  (feet) 

Z — ZO  value  in  the  internal  coordinate  system  (feet) 

V — VO  value  in  the  internal  coordinate  system  (ft/sec) 

GM — GMO  value  in  the  internal  coordinte  system  (radians) 

CHI' — CHIO  value  in  the  internal  coordinate  system  (radians) 

T — TO  value  from  the  command  input  file  (seconds). 

SUBROUTINE  F 

The  subroutine  F  computes  the  instantaneous  time  rates  of  change  of  an  aircraft’s 
Cartesian  position,  speed,  pitch  and  yaw  angles,  and  mass.  Subroutine  F  is  called  by 
subroutine  RK  four  times  to  integrate  these  time  rates  over  the  interval  DT.  The  passed 
parameters  are: 

X — the  X  value  at  one  of  four  integration  steps 
Y — the  Y  value  at  one  of  four  integration  steps 
Z — the  Z  value  at  one  of  four  integration  steps 
V — the  V  value  at  one  of  four  integration  steps 
GM — the  y  value  at  one  of  four  integration  steps 
CHI — the  %  value  at  one  of  four  integration  steps 
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MAC — the  value  of  the  aircraft's  mass,  ACM,  at  one  of  four  integration  steps 
T — the  present  simulated  time 
F's  local  variables  are: 

H — altitude  of  aircraft 

ROH — the  atmospheric  density  p  at  altitude  H 

VSQ — square  of  present  velocity,  V2 

ALNG — longitudinal  acceleration 

APTCH — pitch  acceleration 

AYW — yaw  acceleration 

ATOTAL — magnitude  of  the  total  acceleration. 

FUNCTION  FALPHA 

The  function  FALPHA  is  called  by  subroutine  AERO.  This  function  filters  the  angle 
of  attack.  A  first-order  linear  time  filter  is  used  in  the  calculation.  This  filter  is  discussed  in 
detail  in  [A- 1  ] .  The  passed  parameter  is: 

ALPHAI— commanded  angle  of  attack,  oq. 

FALPHA's  local  variables  are: 

TERM — the  exponential  e'DT/A0ATAEf 

FALPHA — function  value  returned. 

FUNCTION  FCHIDOT 

The  function  FCHIDOT  is  called  from  the  various  control  routines.  This  function 
filters  the  yaw  rate.  A  first-order  linear  time  filter  is  used  in  the  calculation.  This  filter  is 
discussed  in  detail  in  [A-l],  The  passed  parameter  is: 

CHIDOT — commanded  yaw  rate,  X\. 

FCHIDOT's  local  variables  are: 

TERM — the  exponential  e-DELTA/C HITMJ 
FCHIDOT — function  value  returned. 
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FUNCTION  FDF1X 


The  function  FDF1X  is  called  by  subroutine  NWTN1D.  It  computes  the  derivative 
of  the  equation  found  in  function  FUNCF1  with  respect  to  the  independent  variable.  In 

this  case,  the  independent  variable  is  the  instructed  angle  of  attack,  oq.  The  passed 
parameter  is: 


ALPFIAI — the  instructed  angle  of  attack,  ai. 

FDFIX's  local  variable  is: 

FDF1X — function  value  returned. 

FUNCTION  FGMDOT 

The  function  FGMDOT  is  called  from  the  various  control  subroutines.  This  function 
filters  the  pitch  rate.  A  first-order  linear  time  filter  is  used  in  the  calculation.  This  filter  is 
discussed  in  detail  in  [A- 1  ] .  The  passed  parameter  is: 

GMDOTI — the  commanded  pitch  velocity,  Yj, 

FGMDOT’s  local  variables  are: 

TERM — the  exponential  e~DELTA/GMTAU 
FGMDOT — function  value  returned. 

FUNCTION  FMU 

The  function  FMU  is  called  from  subroutine  AERO.  This  function  filters  the  roll 
angle.  A  first-order  linear  time  filter  is  used  in  the  calculation.  This  filter  is  discussed  in 
detail  in  [A- 1  ] .  The  passed  parameter  is: 

MUI — commanded  velocity  roll  angle,  |lp 
FMU's  local  variables  are: 

TERM — the  exponential  e'OT/MUTAU 
FMU — function  value  returned. 

FUNCTION  FUNCF1 

The  function  FUNCF1  is  called  by  the  subroutine  NWTN1D.  This  function 
calculates  the  difference  between  the  commanded  pitch  rate  and  the  computed  pitch  rate 
based  on  the  aircraft's  present  aerodynamic  conditions.  Subroutine  NWTN1D  uses  this 


A- 1 1 


equation  to  estimate  a  value  of  the  instructed  angle  of  attack,  oq,  that  will  minimize  the 
difference.  The  passed  parameter  is: 

ALPHAI — the  passed  X  value  from  NWTN1D,  which  is  the  instructed  angle 
of  attack,  oq 

FUNCFl's  local  variable  is: 

FUNCF1 — function  value  returned. 

FUNCTION  FVDOT 

The  function  FVDOT  is  called  from  the  various  control  routines.  This  function 
filters  the  longitudinal  acceleration.  A  first-order  linear  time  filter  is  used  in  the  calculation. 
This  filter  is  discussed  in  detail  in  [A-l],  The  passed  parameter  is: 


VDOTI — the  commanded  longitudinal  acceleration,  Vi 

FUDOT's  local  variables  are: 

TERM— the  exponential  e'DELTA/VTAU 

FVDOT — function  value  returned. 

SUBROUTINE  NWTN1D 

The  subroutine  NWTN1D  is  called  from  subroutine  AERO.  NWTN1D  performs  a 
one-dimensional  root-finding  algorithm  using  the  Newton-Raphson  method.  This 
algorithm  can  be  found  in  [A-2].  The  passed  parameter  is: 

X — the  passed  parameter  from  subroutine  AERO;  in  this  case,  X  is  equal  to 
the  instructed  angle  of  attack,  oq. 

NWTNID's  local  variables  are: 

FLAG — logical  variable  for  the  DO  WHILE  loop 

EPSLNX — degree  of  tolerance  for  the  difference  between  the  present  value  of 
X  and  the  value  from  the  most  recent  iteration  of  the  root-finding 
algorithm,  XNEW;  once  this  difference  is  less  than  the  value  of  EPSLNX, 
the  DO  WHILE  loop  ends  and  the  value  of  X  is  returned  to  AERO 

JCOUNT — an  integer  counter  keeping  track  of  the  number  of  iterations  that 
the  Newton-Raphson  method  is  taking  to  satisfy  the  tolerance  EPSLNX 

FUNC1 — function  value  returned  from  FUNCF1 

DF1X — function  value  returned  from  FDF1X 
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DET — the  determinant  of  the  matrix  derived  from  the  system  of  equations 
being  used  in  the  root-finding  algorithm;  in  the  one-dimensional  case,  DET 
is  equal  to  DF1X 

DETINV — the  inverse  of  the  determinent,  1/DET 
XNEW — the  computed  value  of  the  root. 

FUNCTION  P 

The  function  P  computes  the  standard  pressure  as  a  function  of  elevation,  H.  The 
passed  parameter  is: 

H — altitude  above  sea  level  (feet). 

P's  local  variable  is: 

P — pressure  at  altitude  H  (lb/ft2). 

SUBROUTINE  PARDER 

The  subroutine  PARDER  is  called  from  program  MOFO.  A  two-dimensional  table 
ARRAY  and  an  independent  variable  array  X  are  received.  Subroutine  PARDER  computes 
partial  derivatives  at  each  position  in  ARRAY  with  respect  to  the  independent  variable 
array  X.  The  method  used  is  the  central-difference  formula,  except  at  the  endpoints.  The 
forward-difference  formula  is  performed  on  the  first  endpoint  while  the  backward- 
difference  formula  is  performed  on  the  last  element  in  the  array.  The  passed  parameters 
are: 


ARRAY(NMAX,NMAX) — the  two-dimensional  table  that  is  passed  to 
PARDER 

DERIV(NMAX,NMAX) — the  returned  two-dimensional  table  of  partial 
derivatives  of  ARRAY (NMAX,NMAX)  with  respect  to  X(NMAX) 

X(NMAX) — the  one-dimensional  independent  variable  array  that  the 
derivative  will  be  taken  with  respect  to 

IMAX — the  maximum  number  of  values  in  the  array  X(NMAX),  also 
corresponds  to  the  number  of  rows  in  ARRAY(NMAX,NMAX) 

JMAX — the  number  of  columns  of  ARRAY(NMAX,NMAX). 

PARDER’s  local  variables  are: 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 

I — row  index  for  the  DO  loop 

J — column  index  for  the  DO  loop. 
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FUNCTION  RHO 


The  function  RHO  computes  the  standard  atmospheric  density  p  as  a  function  of 
elevation,  H.  The  passed  parameter  is: 

H — altitude  above  sea  level  (feet). 

RHO's  local  variable  is: 

RHO — atmospheric  density  as  a  function  of  height  (slugs/ft3). 

SUBROUTINE  RK 

The  subroutine  RK  performs  fourth-order  Runge-Kutta  integration  of  a  system  of 
first-order  nonlinear  differential  equations  (i.e.,  the  equations  of  motion).  The  algorithm 
uses  a  time  step  DT.  Subroutine  F  is  called  four  times  to  complete  the  integration.  See 
[A-l]  and  [A-2]  for  a  more  detailed  explanation  on  Runge-Kutta  integration.  The  passed 
parameters  are: 

X — the  present  X  coordinate  of  the  aircraft  passed  from  a  control  routine;  at 
the  end  of  RK,  X  is  updated  over  the  time  interval  DT 

Y — the  present  Y  coordinate  of  the  aircraft  passed  from  a  control  routine;  at 
the  end  of  RK,  Y  is  updated  over  the  time  interval  DT 

Z — the  present  Z  coordinate  of  the  aircraft  passed  from  a  control  routine;  at  the 
end  of  RK,  Z  is  updated  over  the  time  interval  DT 

V — the  present  magnitude  of  the  velocity  vector  passed  from  a  control  routine; 
at  the  end  of  RK,  V  is  updated  over  the  time  interval  DT 

GM — the  present  pitch  angle  of  the  velocity  vector  passed  from  a  control 
routine;  at  the  end  of  RK,  GM  is  updated  over  the  time  interval  DT 

CHI — the  present  yaw  angle  of  the  velocity  vector  passed  from  a  control 
routine;  at  the  end  of  RK,  CHI  is  updated  over  the  time  interval  DT 

T — the  present  time  of  the  simulation  passed  from  a  control  routine;  at  the  end 
of  RK,  T  is  updated  over  the  time  interval  DT. 

RK’s  local  variables  are: 


FI  1 — X;  the  velocity  component  in  the  X  direction  after  the  first  call  to 
subroutine  F 


F21 — Y;  the  velocity  component  in  the  Y  direction  after  the  first  call  to 
subroutine  F 


F3 1 — Z;  the  velocity  component  in  the  Z  direction  after  the  first  call  to 
subroutine  F 
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F41 — V;  the  acceleration  component  along  the  velocity  vector  after  the  first 
call  to  subroutine  F 


F5 1 — y  the  pitch  rate  after  the  first  call  to  subroutine  F 


F61 — %;  the  yaw  rate  after  the  first  call  to  subroutine  F 


F7 1 — m;  the  rate  of  change  in  the  aircraft's  mass  after  the  first  call  to 
subroutine  F 

XSIN — value  as  a  function  of  X,  DT,  and  FI  1  which  is  passed  to  subroutine 
F  during  the  second  call 

YSIN — value  as  a  function  of  Y,  DT,  and  F21,  which  is  passed  to  subroutine 
F  during  the  second  call 

ZSIN— value  as  a  function  of  Z,  DT,  and  F31,  which  is  passed  to  subroutine 
F  during  the  second  call 

VSIN — value  as  a  function  of  V,  DT,  and  F41,  which  is  passed  to  subroutine 
F  during  the  second  call. 

GMSIN — value  as  a  function  of  GM,  DT,  and  F51,  which  is  passed  to 
subroutine  F  during  the  second  call 

CHSIN — value  as  a  function  of  CHI,DT,  and  F61,  which  is  passed  to 
subroutine  F  during  the  second  call 

ACMSIN — value  as  a  function  of  ACM,  DT,  and  F71,  which  is  passed  to 
subroutine  F  during  the  second  call 


F12 — X;  the  velocity  component  in  the  X  direction  after  the  second  call  to 
subroutine  F 


F22 — Y;  the  velocity  component  in  the  Y  direction  after  the  second  call  to 
subroutine  F 


F32 — Z;  the  velocity  component  in  the  Z  direction  after  the  second  call  to 
subroutine  F 


F42 — V;  the  acceleration  component  along  the  velocity  vector  after  the  second 
call  to  subroutine  F 
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F52 — %  the  pitch  rate  after  the  second  call  to  subroutine  F 


F62 — %;  the  yaw  rate  after  the  second  call  to  subroutine  F 


F72 — m;  the  rate  of  change  in  the  aircraft's  mass  after  the  second  call  to 
subroutine  F 

XSSIN — value  as  a  function  of  X,  DT,  and  F12  which  is  passed  to 
subroutine  F  during  the  third  call 

YSSEN — value  as  a  function  of  Y,  DT,  and  F22,  which  is  passed  to 
subroutine  F  during  the  third  call 

ZSSIN — value  as  a  function  of  Z,  DT,  and  F32,  which  is  passed  to 
subroutine  F  during  the  third  call 

VSSIN — value  as  a  function  of  V,  DT,  and  F42,  which  is  passed  to 
subroutine  F  during  the  third  call 

GMSSIN — value  as  a  function  of  GM,  DT,  and  F52,  which  is  passed  to 
subroutine  F  during  the  third  call 

CHS  SIN — value  as  a  function  of  CHI,  DT,  and  F62,  which  is  passed  to 
subroutine  F  during  the  third  call 

ACMSSIN — value  as  a  function  of  ACM,  DT,  and  F72,  which  is  passed  to 
subroutine  F  during  the  third  call 


F13 — X;  the  velocity  component  in  the  X  direction  after  the  third  call  to 
subroutine  F 


F23 — Y;  the  velocity  component  in  the  Y  direction  after  the  third  call  to 
subroutine  F 


F33 — Z;  the  velocity  component  in  the  Z  direction  after  the  third  call  to 
subroutine  F 


F43 — V;  the  acceleration  component  along  the  velocity  vector  after  the  third 
call  to  subroutine  F 


F53 — y,  the  pitch  rate  after  the  third  call  to  subroutine  F 


F63 — %;  the  yaw  rate  after  the  third  call  to  subroutine  F 
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F73 — m;  the  rate  of  change  in  the  aircraft's  mass  after  the  third  call  to 
subroutine  F 

XSP1 — value  as  a  function  of  X,  DT,  and  F13,  which  is  passed  to  subroutine 
F  during  the  fourth  call 

YSP1 — value  as  a  function  of  Y,  DT,  and  F23,  which  is  passed  to  subroutine 
F  during  the  fourth  call 

ZSP1 — value  as  a  function  of  Z,  DT,  and  F33,  which  is  passed  to  subroutine 
F  during  the  fourth  call 

VSP1 — value  as  a  function  of  V,  DT,  and  F43,  which  is  passed  to  subroutine 
F  during  the  fourth  call 

GMSP1 — value  as  a  function  of  GM,  DT,  and  F53,  which  is  passed  to 
subroutine  F  during  the  fourth  call 

CHSP1 — value  as  a  function  of  CHI,  DT,  and  F63,  which  is  passed  to 
subroutine  F  during  the  fourth  call 

ACMSP1— value  as  a  function  of  ACM,  DT,  and  F73,  which  is  passed  to 
subroutine  F  during  the  fourth  call 


F14 — X;  the  velocity  component  in  the  X  direction  after  the  fourth  call  to 
subroutine  F 


F24 — Y;  the  velocity  component  in  the  Y  direction  after  the  fourth  call  to 
subroutine  F 


F34 — Z;  the  velocity  component  in  the  Z  direction  after  the  fourth  call  to 
subroutine  F 


F44 — V;  the  acceleration  component  along  the  velocity  vector  after  the  fourth 
call  to  subroutine  F 


F54 — y;  the  pitch  rate  after  the  fourth  call  to  subroutine  F 

¥64 — %;  the  yaw  rate  after  the  fourth  call  to  subroutine  F 

¥14 — m;  the  rate  of  change  in  the  aircraft's  mass  after  the  fourth  call  to 
subroutine  F. 
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SUBROUTINE  SPLIE2 

The  subroutine  SPLIE2  is  called  from  program  MOFO.  Given  an  M  by  N  tabulated 
function  YA,  and  tabulated  independent  variables  XI A  (M  values)  and  X2A  (N  values), 
this  routine  constructs  one-dimensional  natural  cubic  splines  of  the  rows  of  YA  and  returns 
the  second  derivatives  in  the  array  Y2A.  The  FORTRAN  code  was  extracted  from  [A-3]. 
The  passed  parameters  are: 

XIA(NMAX) — the  passed  one-dimensional  array  containing  independent 
variables  of  size  M 

X2A(NMAX>— -the  passed  one-dimensional  array  containing  independent 
variables  of  size  N 

YA(NMAX,NMAX) — the  passed  two-dimensional  table  containing  the 
tabulated  function  values  of  XIA(NMAX)  and  X2A(NMAX)  of  size  M  by 
N 

M — integer  describing  the  actual  size  of  XI  A(NMAX) 

N — integer  describing  the  actual  size  of  X2A(NMAX) 

Y2A(NMAX,NMAX) — the  returned  two-dimensional  table  containing  the 
second-derivatives  of  the  tabulated  function  values  in  YA(NMAX,NMAX). 

SPLB^'s  local  variables  are: 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 

YTMP(NMAX) — one-dimensional  array  to  store  all  the  values  of  a  row  of 
Y  A(NMAX,NMAX) 

Y2TMP(NMAX) — the  second  derivatives  for  each  value  in  YTMP(NMAX) 
returned  when  subroutine  SPLINE  is  called 

J — counter  for  DO  loop  and  used  as  an  array  index  term 

K — counter  for  DO  loop  and  used  as  an  array  index  term. 

SUBROUTINE  SPLIN2 

The  subroutine  SPLIN2  is  called  by  various  routines.  Given  XI  A,  X2A,  YA,  M,  N 
as  described  in  SPLIE2  and  Y2A  as  produced  by  that  routine;  and  given  a  desired 
interpolating  point  XI,  X2;  this  routine  returns  an  interpolated  function  value  Y  by  bicubic 
spline  interpolation.  The  FORTRAN  code  was  extracted  from  [A-3].  The  passed 
parameters  are: 

XIA(NMAX) — the  array  as  described  in  subroutine  SPLIE2 
X2A(NMAX) — the  array  as  described  in  subroutine  SPLEE2 
YA(NMAX,NMAX) — the  table  as  described  in  subroutine  SPLIE2 
Y2A(NMAX,NMAX) — the  table  as  described  in  subroutine  SPLIE2 


A-18 


M — integer  as  described  in  subroutine  SPLIE2 
N — integer  as  described  in  subroutine  SPLIE2 
XI — desired  interpolating  point  for  XIA(NMAX) 

X2 — desired  interpolating  point  for  X2A(NMAX) 

Y — the  returned  interpolated  function  value. 

SPLIN2's  local  variables  are: 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 

YTMP(NMAX)— one-dimensional  array  to  store  all  the  values  of  a  row  of 
Y  A(NMAX,NM  AX) 

Y2TMP(NMAX) — one-dimensional  array  to  store  all  the  values  of  a  row  of 
Y2A (NMAX, NMAX),  then  used  to  store  the  second  derivatives  of  the 
array  YYTMP(NMAX)  when  SPLINE  is  called 

YYTMP(NMAX) — one-dimensional  array  created  by  interpolating 
Y A(NMAX,NMAX)  at  X2  for  all  the  rows 

J— -counter  for  DO  loop  and  used  as  an  array  index  term 

K — counter  for  DO  loop  and  used  as  an  array  index  term. 

SUBROUTINE  SPLINE 

The  subroutine  SPLINE  is  called  by  program  MOFO.  Given  arrays  X  and  Y  of 
length  N  containing  a  tabulated  function,  and  given  value  YP1  and  YPN  for  the  first 
derivative  of  the  interpolating  function  at  points  1  and  N,  respectively,  this  routine  returns 
an  array  Y2  of  length  N,  which  contains  the  second  derivatives  of  the  interpolating  function 
at  the  tabulated  points  Xj.  If  YP1  and/or  YPN  are  equal  to  1  x  1030  or  larger,  the  routine  is 
signalled  to  set  the  corresponding  boundary  condition  for  a  natural  spline,  with  zero  second 
derivative  on  that  boundary.  The  FORTRAN  code  was  extracted  from  [A-3].  The  passed 
parameters  are: 

X(NMAX) — the  passed  one-dimensional  array  containing  the  tabulated 
independent  variables  of  size  N 

Y(NMAX) — the  passed  one-dimensional  array  containing  the  tabulated 
independent  function  values  of  size  N 

N — integer  describing  the  actual  size  of  X(NMAX)  and  Y  (NMAX) 

YP1 — first  derivative  of  the  interpolating  function  at  point  1 

YPN — first  derivative  of  the  interpolating  function  at  point  N 


A- 19 


Y2(NMAX) — one-dimensional  array  containing  the  second  derivatives  of  the 
interpolating  function  at  the  tabulated  points  Xj. 

SPLINE'S  local  variables  are: 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 

U(NMAX) — an  array  that  acts  as  an  intermediate  step  in  the  calculation  of  the 
second  derivative 

I — counter  for  DO  loop  and  used  as  an  array  index  term 
K — counter  for  DO  loop  and  used  as  an  array  index  term. 

SUBROUTINE  SPLINT 

The  subroutine  SPLINT  is  called  by  various  subroutines.  Given  the  arrays  XA  and 
YA  of  length  N,  which  tabulate  a  function  (with  the  XAjS  in  order),  given  the  array  Y2A, 
which  is  the  output  from  SPLINE  above,  and  given  a  value  of  X,  this  routine  returns  a 
cubic-spline  interpolated  value  Y.  The  FORTRAN  code  was  extracted  from  [A-3].  The 
passed  parameters  are: 

XA(NMAX) — the  array  as  described  in  subroutine  SPLINE 
YA(NMAX) — the  array  as  described  in  subroutine  SPLINE 
Y2A(NMAX — the  array  as  described  in  subroutine  SPLINE 
N — integer  as  described  in  subroutine  SPLINE 
X — desired  interpolating  point  of  the  array  XA(NMAX) 

Y — the  returned  interpolated  function  value. 

SPLINT's  local  variables  are: 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 

KLO — index  value  to  find  the  largest  value  of  XA(NMAX)  that  still  satisfies 
the  relationship,  XA(KLO)  <  X 

KHI — index  value  to  find  the  smallest  value  of  XA(NMAX)  that  still  satisfies 
the  relationship,  XA(KHI)  >  X 

K — index  value  to  locate  KHI  and  KLO 

H — difference  between  XA(KHI)  and  XA(KLO) 

A— the  ratio:  [XA(KHI)  -  X]/H 

B— the  ratio:  [X  -  XA(KLO)]/H. 
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SUBROUTINE  TABIN 


The  subroutine  TABIN  is  called  from  program  MOFO  to  read  in  the  power  and 
aerodynamic  tables  from  the  aircraft  characteristic  file.  A  flag  variable  is  passed  to  indicate 
whether  the  file  is  composed  of  real  or  integer  values.  Integer  inputs  are  converted  to  real. 
Variable  SCALE  multiplies  the  table  values,  if  not  equal  to  one.  The  passed  parameters  are: 

XIARGS(NMAX) — the  independent  array  values  associated  with  each  of  the 
rows  of  the  two-dimensional  table 

X2ARGS(NMAX) — the  independent  array  values  associated  with  each  of  the 
columns  of  the  two-dimensional  table 

TABLE(NMAX,NMAX) — table  to  be  read  in 

NROWS — number  of  rows  in  table 

NCOLS — number  of  columns  in  table 

KFLAG — flag  that  indicates  whether  the  table  and  the  first  argument  is 
composed  of  integer  or  real  values. 

TABIN's  local  variables  are: 

NMAX — parameter  descriptive  of  the  maximum  array  size  currently  set  to  14 
IARGS(NMAX)— integer  form  of  the  X 1 ARGS (NMAX)  array 
ITABLE(NMAX,NMAX) — integer  table  to  be  read  in 
SCALE — a  scale  factor  that  multiplies  the  table 
IR — counter  for  number  of  rows  in  the  DO  loop 
IC — counter  for  number  of  columns  in  the  DO  loop. 

SUBROUTINE  WAYPT 

The  subroutine  WAYPT  is  called  by  subroutine  CONTROL  to  construct  a  flight 
segment.  The  present  position  and  velocity  vector  of  the  aircraft  are  passed  in  this  call 
statement.  The  subroutine  then  reads  in  its  commanded  position  and  calculates  the 
directional  change  necessary  to  point  the  velocity  vector  to  the  desired  waypoint.  Finally,  it 
constructs  a  flight  path  by  calling  subroutine  RK  to  achieve  the  desired  position.  The 
passed  parameters  are: 

X — the  present  X  coordinate  of  the  aircraft 

Y — the  present  Y  coordinate  of  the  aircraft 

Z — the  present  Z  coordinate  of  the  aircraft 

V — the  present  magnitude  of  the  velocity  of  the  aircraft 

GM — the  present  velocity  pitch  angle 
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CHI — the  present  velocity  yaw  angle 
T — the  present  simulated  time  of  the  model. 

WAYPT's  local  variables  are: 

FLAG — logical  for  the  DO  loop 
XF — X  value  of  desired  waypoint 
ZF — Z  value  of  desired  waypoint 
VF — V  value  of  desired  waypoint 

DELX — difference  between  present  X  value  and  desired  X  value 
DELY — difference  between  present  Y  value  and  desired  Y  value 
DELZ — difference  between  present  Z  value  and  desired  Z  value 
VX — the  component  of  the  velocity  vector  in  the  X  direction 
VY — the  component  of  the  velocity  vector  in  the  Y  direction 
VZ — the  component  of  the  velocity  vector  in  the  Z  direction 
RANGE — slant  range  between  present  position  and  desired  waypoint 

OMGALOSX — ^losx;  see  [A-l]  for  details 

OMGALOSY — ^LOSyi see  IA-1]  for  details 

OMGALOSZ — £2l0Sz;  see  [A-l]  for  details 

DELV — the  difference  between  VF  and  V 

AVEV — average  speed  of  VF  and  V 

TTI — time  to  intercept  calculated  by  RANGE  and  AVEV 

DRAGACCEL — maximum  acceleration  due  to  drag 

RATIO — variable  used  as  a  check  to  be  sure  that  the  commanded  directional 
changes  are  not  greater  than  the  user-imposed  g-limits;  if  they  are,  they  are 
scaled  accordingly 

CHIDOTI — the  computed  yaw  rate  command  before  filtering, 
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GMDOTI — the  computed  pitch  rate  command  before  filtering,  y\ 

VDOTI — the  computed  longitudinal  acceleration  command  before  filtering,  Vi 
TF — simulated  time  when  the  acceleration  commands  are  updated. 


A-23 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REFERENCES 


[A-l]  CNA  Research  Memorandum  89-63,  An  Analyst's  Guide  to  the  MOFO  Model ,  by 
Michael  W.  Price  and  James  E.  Schliessmann,  Unclassified,  May  1989 

(27 890063 11 

[A-2]  Robert  W.  Hombeck.  Numerical  Methods.  New  Jersey:  Quantum  Publishers, 
1975 

[A- 3]  William  H.  Press  et.  al.  Numerical  Recipes.  New  York:  Cambridge  University 
Press,  1986 


1.  The  number  in  parentheses  is  a  CNA  internal  control  number. 


APPENDIX  B 
SOURCE  CODE 


VDOO-JO'iUl.C.UJtO  U  to 


c 


c 


PARAMETER  (NMAX=14) 

IMPLICIT  REAL  *  8 { A— H , 0— Z ) 

REAL*8  KAY , MACH , MTOFT , MU , MUO , MUTAU , MV , MG 

REAL*8  MACHCD , MACHCL , MCHAB , MCHMI L 

COMMON  /AERO/A I RS PD ( NMAX ) , AOAMAX ( NMAX ) , AAMX S D ( NMAX ) , 

2  CLFTCD ( NMAX ), MACHCD (NMAX) , 

3  CD, CDRAG( NMAX, NMAX) , CDRGSD ( NMAX , NMAX ) , 

4  AOACL ( NMAX ), MACHCL (NMAX) , 

5  CL , CL I FT ( NMAX, NMAX ) , CLFTSD ( NMAX , NMAX ) , 

6  DERCD,DCDDCL( NMAX, NMAX) , DCDDCLSD( NMAX , NMAX ) , 

7  DERCL,DCLDA( NMAX, NMAX) , DC LDAS D ( NMAX , NMAX ) 

COMMON  /AIRCFT/ACM, AREA, BINGO, WTAC , WTFLMX , WTFUEL 

COMMON  /AOA/AL P HA , AL  P HA 0 , AOATAU 

COMMON  /CONSTS/G , COUNT , I COUNT , MTOFT , PO , PI , RHOO , ROOT , RTODEG 

COMMON  /CONTRL/CHIDOT, CHIDOTO , CHITAU , 

GMDOT , GMDOTO , GMTAU , 

VDOT , VDOT0 , VTAU , 

DELTA, DT, GEES, KAY 

COMMON  /CYCLE/MV, MG, REFOR, SNG, CSG, SNCHI ,CSCHI , CT 

COMMON  /FNCTN/F1 , F2 , F3 , F4 , F5 , F6 , F7 

COMMON  /POWER/ALT (NMAX) , SPDMAX ( NMAX ) , SPMXSD ( NMAX ) , 
ALTMIL(NMAX) , MCHMI L ( NMAX ) , 

PWRMIL( NMAX, NMAX) , PRMLSD ( NMAX , NMAX ) , 

FULMIL( NMAX, NMAX) , FLMLSD ( NMAX , NMAX ) , 

ALTAB ( NMAX ) , MCHAB { NMAX ) , 

PWRAB ( NMAX , NMAX ) , PRABSD ( NMAX , NMAX ) , 

FULAB ( NMAX , NMAX ) , FLABSD ( NMAX , NMAX ) , 

PCTRPM ( NMAX ) ,TSFC(NMAX) , TSFCSD ( NMAX ) , 

MACH , TH , THMAX 

COMMON  /ROLL/MU , MUO , MUTAU 

COMMON  /TBLMAX/NMCHMX , NAOAMX , NALTML , NMCHML , NALTAB , NMCHAB , 

2  NCLCD , NMCHCD , NAOACL , NMCHCL , NCTRPM 
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PROGRAM  MOFO 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  PROGRAM:  MODEL  OF  FLYING  OBJECTS  (MOFO)  C 
C  C 
C  AUTHORS:  Michael  W.  Price  C 
C  James  E.  Schliessmann  C 
C  C 
C  CENTER  FOR  NAVAL  ANALYSES  C 
C  4401  FORD  AVENUE  C 
C  P.O.  BOX  16268  C 
C  ALEXANDRIA,  VA  22302-0268  C 
C  (703)  824-2693  C 
C  C 
C  PROGRAM  MOFO  IS  THE  MAIN  ROUTINE  FOR  THE  MOFO  MODEL.  IT  FIRST  C 
C  INITIALIZES  SEVERAL  OF  THE  VARIABLES  USED  THROUGHOUT  THE  REST  OF  THE  C 
C  PROGRAM.  IT  THEN  CONTINUES  ON  TO  READ  IN  THE  AIRCRAFT  DATA  FILE.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  ' CMNBLK/LIST' 

C 

INTEGER  I ALT ( NMAX ) 

C 

DATA  CHIDOTO/O .DO/ 

DATA  G/32.17D0/ 

DATA  GMDOTO/O . DO/ 

DATA  MTOFT/3 . 28084 0D0/ 

DATA  P0/2116 . 2D0/ 

DATA  PI/3. 141 5926535897 9D0/ 

DATA  RHO0/2 . 377D-3/ 

DATA  VDOT0/0.D0/ 

C 

ROOT  =  1.DO/3.5DO 
RTODEG  =  180. DO/PI 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  C 

C  BEGIN  READING  THE  AIRCRAFT  DATA  FILE.  START  BY  READING  IN  THE  C 
C  ARRAY  OF  MAXIMUM  MACH  NUMBER  AS  A  FUNCTION  OF  ALTITUDE  AND  THE  ARRAY  C 
C  OF  MAXIMUM  ANGLE  OF  ATTACK  AS  A  FUNCTION  OF  MACH.  FOR  EACH  ARRAY,  C 
C  SUBROUTINE  SPLINE  IS  CALLED  TO  CREATE  AN  ARRAY  OF  SECOND  DERIVATIVES  C 
C  TO  BE  USED  LATER  WHEN  INTERPOLATING  FOR  SPECIFIC  VALUES.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

READ (1,1001)  NMCHMX , NAOAMX 
C 

READ( 1,1003)  (IALT(I) ,1=1, NMCHMX) 

READ (1,1002)  (SPDMAX(I) ,1=1, NMCHMX) 

DO  69  I  =  1, NMCHMX 

ALT(I)  =  DFLOAT( IALT( I ) ) 

69  CONTINUE 

CALL  SPLINE ( ALT, SPDMAX, NMCHMX, 1 . D3 0 , 1 . D3 0 , SPMXSD ) 

C 

READ (1,1002)  (AIRSPD(I) ,1=1, NAOAMX) 

READ( 1,1002)  (AOAMAX( I) ,1=1, NAOAMX) 

DO  9696  1=1, NAOAMX 

AOAMAX ( I ) =AOAMAX ( I ) /RTODEG 
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9696  CONTINUE 

CALL  SPLINE ( AIRSPD , AOAMAX , NAOAMX, l.D30,l.D30 , AAMXSD ) 

C 

WRITE ( 7 , 3333 ) 

WRITE (8,3334) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  c 
C  READ  IN  THE  2-DIMENSIONAL  ARRAY  OF  MILITARY  POWER  AS  A  FUNCTION  C 
C  OF  ALTITUDE  AND  MACH  NUMBER  BY  CALLING  SUBROUTINE  TAB I N .  CALL  C 
C  SUBROUTINE  SPLIE2  TO  CREATE  AN  ARRAY  OF  SECOND  DERIVATIVES  NEEDED  C 
C  IN  THE  BICUBIC  SPLINE  INTERPOLATION  ROUTINE.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  TAB I N ( ALTM I L , MCHMIL , PWRMIL , NALTML , NMCHML , 2 ) 

CALL  SPLIE2 ( ALTMI L , MCHMIL , PWRMIL , NALTML , NMCHML , PRMLSD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC.CCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  IN  THE  2-DIMENSIONAL  ARRAY  OF  AFTERBURNER  POWER  AS  A  C 
C  FUNCTION  OF  ALTITUDE  AND  MACH  NUMBER.  CREATE  AN  ARRAY  OF  ITS  SECOND  C 
C  DERIVATIVES.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  TAB  I  N't  ALTAB  ,  MCHAB ,  PWRAB  ,  NALTAB  ,  NMCHAB  ,  2  ) 

CALL  SPLIE2 ( ALTAB , MCHAB , PWRAB , NALTAB , NMCHAB , PRABSD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  IN  THE  2-DIMENSIONAL  ARRAY  OF  FUEL  CONSUMPTION  UNDER  C 
C  MILITARY  POWER  AS  A  FUNCTION  OF  ALTITUDE  AND  MACH  NUMBER.  CREATE  C 
C  AN  ARRAY  OF  ITS  SECOND  DERIVATIVES.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  TAB IN ( ALTMI L , MCHMIL , FULMI L , NALTML , NMCHML , 2 ) 

CALL  SPLI E2 ( ALTMI L , MCHMI L , FULMIL , NALTML , NMCHML , FLMLSD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  IN  THE  2-DIMENSIONAL  ARRAY  OF  FUEL  CONSUMPTION  UNDER  C 
C  AFTERBURNER  POWER  AS  A  FUNCTION  OF  ALTITUDE  AND  MACH  NUMBER.  CREATE  C 
C  AN  ARRAY  OF  ITS  SECOND  DERIVATIVES.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  TABIN( ALTAB , MCHAB , FULAB , NALTAB , NMCHAB , 2 ) 

CALL  SPLIE2 ( ALTAB , MCHAB , FULAB , NALTAB , NMCHAB , FLABSD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  IN  THE  2-DIMENSIONAL  ARRAY  OF  DRAG  COEFFICIENT  AS  A  C 
C  FUNCTION  OF  LIFT  COEFFICIENT  AND  MACH  NUMBER.  CREATE  AN  ARRAY  OF  C 
C  ITS  SECOND  DERIVATIVES.  CALL  SUBROUTINE  PARDER  TO  CREATE  AN  ARRAY  C 
C  OF  PARTIAL  DERIVATIVES  WITH  RESPECT  TO  THE  LIFT  COEFFICIENT  FOR  C 
C  EACH  MACH  NUMBER.  CREATE  AN  ARRAY  OF  ITS  SECOND  DERIVATIVES.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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CALL  TAB IN ( CLFTCD , MACHCD , CDRAG , NCLCD , NMCHCD , 1 ) 

CALL  SPLIE2 ( CLFTCD , MACHCD , CDRAG, NCLCD , NMCHCD , CDRGSD ) 

CALL  PARDER( CDRAG, CLFTCD, NCLCD, NMCHCD, DCDDCL) 

CALL  SPLIE2 { CLFTCD , MACHCD , DCDDCL , NCLCD , NMCHCD , DCDDCLSD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  c 
C  READ  IN  THE  2-DIMENSIONAL  ARRAY  OF  LIFT  COEFFICIENT  AS  A  C 
C  FUNCTION  OF  ANGLE  OF  ATTACK  AND  MACH  NUMBER.  CREATE  AN  ARRAY  OF  ITS  C 
C  SECOND  DERIVATIVES.  CREATE  AN  ARRAY  OF  ITS  PARTIAL  DERIVATIVES  WITH  C 
C  RESPECT  TO  ANGLE  OF  ATTACK.  THEN  CREATE  AN  ARRAY  OF  ITS  SECOND  C 
C  DERIVATIVES.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  TABIN ( AOACL , MACHCL , CLI FT , NAOACL , NMCHCL , 1 ) 

DO  6969  1=1, NAOACL 

AOACL ( I ) = AOACL ( I ) /RTODEG 
6969  CONTINUE 

CALL  SPLI E2 ( AOACL , MACHCL , CLI FT , NAOACL , NMCHCL , CLFTSD ) 

CALL  PARDER (CLIFT, AOACL , NAOACL , NMCHCL , DCLDA ) 

CALL  SPLI E2 ( AOACL , MACHCL , DCLDA , NAOACL , NMCHCL , DCLDASD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  IN  THE  ONE  DIMENSIONAL  TABLE  OF  THRUST  SPECIFIC  FUEL  C 

C  CONSUMPTION  AS  A  FUNCTION  OF  THE  PERCENT  OF  RPM'S.  C 

C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
READ (1,998)  NCTRPM 

READ( 1,1002)  (PCTRPM(I) ,1=1, NCTRPM) 

READ( 1,1002)  (TSFC(I) ,1=1, NCTRPM) 

CALL  SPLINE ( PCTRPM , TSFC , NCTRPM , 1 . D30 , 1 . D3  0 , TSFCSD ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  C 

c  READ  IN  THE  DIFFERENTIAL  TIME  STEP  SIZE,  THE  DELTA  TIME  STEP  C 
C  SIZE,  AND  THE  NAVIGATION  CONSTANT  FOR  PROPORTIONAL  NAVIGATION.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

READ( 1,999)  DT, DELTA, KAY 
COUNT  =  1.D0/DT 
ICOUNT  =  NINT( COUNT) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  DIFFERENTIAL  TIME  STEP  SIZE  AND  AERODYNAMIC  TIME  CONSTANTS  C 
C  FOR  ANGLE  OF  ATTACK  AND  ROLL  ANGLE  FILTERING.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

READ (1,1000)  AOATAU , MUTAU , GMTAU , CH I TAU , VTAU 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 

C  READ  IN  THE  AIRCRAFT'S  PLANFORM  AREA,  BASE  WEIGHT,  AND  MAXIMUM  C 

C  FUEL  CAPACITY.  THEN  CALL  SUBROUTINE  CONTROL  C 
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c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

READ (1,2000)  AREA , WTAC , WTFLMX 
C 

CALL  CONTROL 
C 

STOP 

C  *** 

998  FORMAT ( I 1 0 ) 

999  FORMAT ( /3 Fl 0 . 0 ) 

1000  FORMAT ( /5 Fl 0 . 0 ) 

1001  FORMAT ( /2 I 1 0 ) 

1002  FORMAT( 10X, 14F5 . 0 ) 

1003  FORMAT( 10X, 1415  ) 

2000  FORMAT(/3F15 . 0 ) 

3333  FORMAT ( IX,  't  (s),  X  (m),  Y  (m),  Z  (m),  V  (m/s),  BODY  YAW  (deg.) 

1,  BODY  PITCH  (deg.),  BODY  ROLL  (deg.)') 

3334  FORMAT  ( IX,  't  (s),  m  (slugs),  GAMMA  (deg.),  CHI  (deg..), 

&  AOA  (deg.),  MU  (deg.),  T  (lbs.),  Al  (m/s**2), 

&  Ap  (m/s**2),  Ay  (m/s**2),  At  (m/s**2)') 

C  *** 

END 
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SUBROUTINE  SPLINE (X,Y,N,YPl,YPN,Y2) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  GIVEN  ARRAYS  X  AND  Y  OF  LENGTH  N  CONTAINING  A  TABULATED  FUNCTION, 

C  AND  GIVEN  VALUES  YPl  AND  YPN  FOR  THE  FIRST  DERIVATIVE  OF  THE 

C  INTERPOLATING  FUNCTION  AT  POINTS  1  AND  N,  RESPECTIVELY,  THIS  ROUTINE 

C  RETURNS  AN  ARRAY  Y2  OF  LENGTH  N  WHICH  CONTAINS  THE  SECOND  DERIVATIVES 

C  OF  THE  INTERPOLATING  FUNCTION  AT  THE  TABULATED  POINTS  Xi .  IF  YPl 

C  AND/OR  YPN  ARE  EQUAL  TO  1.E30  OR  LARGER,  THE  ROUTINE  IS  SIGNALLED  TO 

C  SET  THE  CORRESPONDING  BOUNDARY  CONDITION  FOR  A  NATURAL  SPLINE,  WITH 

C  ZERO  SECOND  DERIVATIVE  ON  THAT  BOUNDARY. 

C 

C  REFERENCE:  NUMERICAL  RECIPES,  W.H.  PRESS,  et  al. 

C  CAMBRIDGE  U.  PRESS, 1986 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IMPLICIT  REAL *8  (A-H,0-Z) 

C 

PARAMETER  ( NMAX= 1 4 ) 

C 

REAL*8  X(NMAX) , Y ( NMAX ) ,Y2(NMAX) ,U(NMAX) 

C 

IF  (YPl  .GT.  . 99D30 )  THEN 
Y2(l)  =  0 . DO 
U ( 1 )  -  0 .DO 
ELSE 

Y2 ( 1 )  =  -0.SD0 

U ( 1 )  =  ( 3 . D0/( X( 2 ) -X( 1 )))*(( Y( 2 ) -Y (1) )/(X(2)-X(l) ) -YPl ) 

END  IF 
C 

DO  11  I  =  2 , N-l 

SIG  =  (X(I)-X(I-l) )/<X(I+l)-X(I~l) ) 

P  =  SIG*Y2(I-1)+2.D0 
Y2  (  I )  =  ( SIG-1 . DO ) /P 

U(I)  =  ( 6 . DO  * ( (Y( I  +  l)-Y( I) ) / ( X ( 1+1) -X( I))-(Y(I)-Y(I-1)) 

*  /(X(I)-X(I-l) ) )/ ( X ( I+1)-X(I— 1) )— SIG*U( 1-1) )/P 

11  CONTINUE 
C 

IF  (YPN  .GT.  . 99D30 )  THEN 
QN  =  0 . DO 
UN  =  0 . DO 
ELSE 

QN  =  0.5D0 

UN  =  ( 3.D0/(X(N)-X(N-1 ) ) ) *(YPN-< Y(N)-Y(N-1 ) ) /( X ( N ) -X ( N-l ) ) ) 

END  IF 
C 

Y2(N)  =  ( UN— QN*U ( N— 1 ) )/( QN* Y2 ( N-l ) +1 . DO ) 

C 

DO  12  K  =  N-l, 1,-1 

Y2(K)  =  Y2(K)*Y2(K+1)+U(K) 

12  CONTINUE 
C 

RETURN 

END 


SUBROUTINE  SPLINT ( XA , YA, Y2A , N , X , Y ) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  GIVEN  THE  ARRAYS  XA  AND  YA  OF  LENGTH  N,  WHICH  TABULATE  A  FUNCTION  C 
C  (WITH  THE  XAi'S  IN  ORDER),  AND  GIVEN  THE  ARRAY  Y2A ,  WHICH  IS  C 
C  THE  OUTPUT  FROM  SPLINE  ABOVE,  AND  GIVEN  A  VALUE  OF  X,  THIS  ROUTINE  C 
C  RETURNS  A  CUBIC-SPLINE  INTERPOLATED  VALUE  Y.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IMPLICIT  REAL *8  (A-H,0-Z) 

C 

PARAMETER  ( NMAX= 1 4 ) 

C 

REAL *8  XA(NMAX) ,YA(NMAX) , Y2A(NMAX) 

C 

KLO  =  1 

KHI  =  N 
C 

DO  WHILE  ( KHI-KLO  .GT.  1) 

K  =  ( KHI+KLO ) /2 
IF  ( XA ( K )  .GT.  X)  THEN 
KHI  =  K 
ELSE 

KLO  »  K 
END  IF 

END  DO 
C 

H  =  XA(KHI)  -  XA(KLO) 

IF  (H  .EQ.  0 . DO )  THEN 
PAUSE  'Bad  XA  input.' 

ENDIF 

C 

A  =  (XA( KHI )-X)/H 

B  =  (X-XA(KLO) )/H 

Y  =  A*YA( KLO ) +B*YA( KHI ) +■ 

*  ( ( A**3-A ) * Y2A( KLO ) + ( B**3~B ) * Y2A( KHI ) ) * ( H* *2 )/6 . DO 

C 

RETURN 

END 


SUBROUTINE  SPLIE2 ( XlA, X2A , YA, M , N, Y2A ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

c  GIVEN  AN  M  BY  N  TABULATED  FUNCTION  YA,  AND  TABULATED  INDEPENDENT  C 
C  VARIABLES  XlA  (M  VALUES)  AND  X2A  (N  VALUES),  THIS  ROUTINE  CONSTRUCTS  C 
C  ONE-DIMENSIONAL  NATRUAL  CUBIC  SPLINES  OF  THE  ROWS  OF  YA  AND  RETURNS  C 
C  THE  SECOND-DERIVATIVES  IN  THE  ARRAY  Y2A.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IMPLICIT  REAL *8  (A-H,0-Z) 

C 

PARAMETER  (NMAX=14) 

C 

REAL* 8  XlA(NMAX) ,X2A(NMAX) , YA ( NMAX , NMAX ) , Y2A( NMAX , NMAX ) , 

*  YTMP ( NMAX ) ,Y2TMP( NMAX) 

C 

DO  13  J  =  1 ,  M 
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DO  11  K  =  1,N 

YTMP(K)  -  YA(J,K) 

11  CONTINUE 

CALL  SPLINE ( X2A, YTMP ,N,1.D30,1.D30, Y2TMP ) 
DO  12  K  =  1 , N 

Y2A(J,K)  =  Y2TMP ( K ) 

12  CONTINUE 

13  CONTINUE 

RETURN 

END 


SUBROUTINE  SPLIN2 ( XlA, X2A, YA , Y2A , M , N , Xl , X2 , Y ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

GIVEN  XlA ,  X2A,  YA,  M,  N  AS  DESCRIBED  IN  SPLIE2  AND  Y2A  AS 
PRODUCED  BY  THAT  ROUTINE;  AND  GIVEN  A  DESIRED  INTERPOLATING  POINT  Xl , 

X2;  THIS  ROUTINE  RETURNS  AN  INTERPOLATED  FUNCTION  VALUE  Y  BY  BICUBIC 
SPLINE  INTERPOLATION. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IMPLICIT  REAL*8  (A-H,0-Z) 

C 

PARAMETER  (NMAX=14) 

C 

REAL*8  XlA(NMAX) ,X2A(NMAX) , YA ( NMAX , NMAX ) , Y2A( NMAX , NMAX ) , 

*  YTMP { NMAX ) , Y2TMP ( NMAX ) , YYTMP ( NMAX ) 

C 

DO  12  J  =  1 , M 
DO  11  K  =  1 , N 

YTMP(K)  =  YA ( J , K ) 

Y2TMP ( K )  =  Y2A( J , K ) 

11  CONTINUE 

CALL  SPLINT ( X2A, YTMP , Y2TMP , N , X2 , YYTMP ( J ) ) 

12  CONTINUE 
C 

CALL  SPLINE (XlA, YYTMP, M, 1 .D30,l ,D30,Y2TMP) 

CALL  SPLINT ( XlA , YYTMP , Y2TMP , M , Xl , Y ) 

C 

RETURN 

END 
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n  o  o  o  o  o  o 


SUBROUTINE  TAB IN ( XlARGS , X2ARGS , TABLE , NROWS , NCOLS , KFLAG ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  TAB IN  IS  CALLED  FROM  PROGRAM  MOFO  TO  READ  IN  THE 
2-DIMENSIONAL  TABLES  FROM  THE  AIRCRAFT  DATA  FILE.  DEPENDING  ON  THE 
VALUE  OF  KFLAG,  TAB IN  WILL  READ  IN  EITHER  INTEGER  OR  REAL  VALUES. 

IF  VARIABLE  SCALE  IS  NOT  EQUAL  TO  ONE,  THEN  THE  ARRAY  IS  MULTIPLIED 
BY  THIS  FACTOR. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

PARAMETER  (NMAX=14) 

C 

IMPLICIT  REAL *8 ( A— H , 0— Z ) 

REAL *8  TABLE (NMAX,NMAX) , XlARGS ( NMAX) , X2ARGS ( NMAX ) 

C 

INTEGER  I TABLE (NMAX, NMAX ) ,IARGS( NMAX) 

C 

READ (1,10 00)  NROWS, NCOLS, SCALE 
READ( 1,1001)  (X2ARGS(IC),  IC=1, NCOLS) 

IF  (KFLAG  .EQ.  1)  THEN 
DO  10  IR  =  1, NROWS 

READ( 1,1002 ) XlARGS ( IR) , ( TABLE ( IR , IC ) ,IC=1, NCOLS) 

10  CONTINUE 

ELSE  IF  (KFLAG  . EQ .  2)  THEN 
DO  20  IR  =  1, NROWS 

READ ( 1,1003 )IARGS( IR) , ( ITABLE( IR , IC ) ,IC=1, NCOLS) 

XlARGS ( IR )  =  DFLOAT ( IARGS ( IR ) ) 

DO  30  IC  =  1, NCOLS 

TABLE ( IR , I C )  =  DFLOAT ( ITABLE( IR, IC) ) 

30  CONTINUE 

20  CONTINUE 

END  IF 

IF  (SCALE  .NE.  1.D0)  THEN 
DO  40  IR-1, NROWS 
DO  50  IC=1, NCOLS 

TABLE ( IR , IC ) =SCALE *TABLE ( IR, IC ) 

50  CONTINUE 

40  CONTINUE 

ENDIF 
RETURN 

C  *** 

1000  FORMAT ( 10X,2I5,30X,F10.0/) 

1001  FORMAT ( 1 OX , 1 4  F  5 . 0 ) 

1002  FORMAT (4X,F5.0,1X,14F5.0) 

1003  FORMAT (4X,I5,1X,14I5) 

C  *** 

END 
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o  o  o  o  o 


SUBROUTINE  PARDER ( ARRAY , X , IMAX , JMAX , DERIV ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  PARDER  IS  CALLED  FROM  PROGRAM  MOFO.  IT  TAKES  THE 
VALUES  IN  ARRAY  AND  CALCULATES  THE  PARTIAL  DERIVATIVES  WITH  RESPECT 
TO  VARIABLE  ARRAY  X. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

PARAMETER  (NMAX=14) 

C 

IMPLICIT  REAL*8(A-H,0-Z ) 

REAL *8  ARRAY ( NMAX , NMAX ) , DERIV( NMAX , NMAX) ,X(NMAX) 

C 

DO  500  J=1 , JMAX 

DERIV ( 1 , J )  =  ( ARRAY ( 2 , J ) -ARRAY { 1 ,  J  ) )/( X( 2 ) -X( 1 )  ) 

DERIV( IMAX, J)  =  ( ARRAY ( IMAX , J ) -ARRAY ( IMAX-1 , J ) ) / 

&  (X( IMAX)-X( IMAX-1 ) ) 

DO  300  1=2, IMAX-1 

DERI V ( I , J )  =  ( ARRAY ( I +1 , J ) —ARRAY (I— 1,J) )/(X(I+l) —X ( I— 1 ) ) 

300  CONTINUE 
500  CONTINUE 
RETURN 
END 
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SUBROUTINE  CONTROL 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  SUBROUTINE  CONTROL  IS  CALLED  FROM  PROGRAM  MOFO.  ITS  TASK  IS  TO  C 
C  PROCESS  THE  FILE  THAT  COMMANDS  THE  AIRCRAFT  TO  PERFORM  TO  THE  USER'S  C 
C  SPECIFICATIONS.  ORDERS  ARE  READ  IN  AND  THEN  TRANSFERRED  TO  THE  C 
C  APPRORIATE  SUBROUTINE  TO  COMPLETE  THE  COMMAND.  AT  THE  END  OF  THE  C 
C  COMMAND  FILE,  THE  SUBROUTINE  KICKS  OUT  THUS  ENDING  A  MOFO  SESSION.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  ' CMNBLK/LIST' 

C 

REAL* 8  MXMACH 
C 

CHARACTER* 2  ORDERS 
C 

LOGICAL  FLAG 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 
C  READ  IN  THE  INITIAL  CONDITIONS.  C 
C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


READ ( 2 , *  ) 
READ ( 2 , * ) 
READ (2,5) 
READ (2,5 ) 
READ ( 2 , 5 ) 
READ ( 2 , * ) 


XO , YO , 20 , VO 

CH10 , GMO , MUO , ALPHAO 

TO , WTFUEL , BINGO , WTARM 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  c 
c  A  CONVERSION  TO  MIELE'S  COORDINATE  SYSTEM  IS  PERFORMED.  A  CHECK  C 
C  ON  FUEL  CAPICITY  IS  MADE.  THE  WEIGHTS  OF  THE  VEHICLE,  ITS  FUEL,  AND  C 
C  EXTERNAL  STORES  ARE  ADDED , AND  CONVERTED  TO  MASS.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

X  =  XO  *MTOFT 

Y  =  -YO*MTOFT 
Z  =  — ZO*MTOFT 

c 

V  =  VO*MTOFT 

GM  =  GMO/RTODEG 

CHI  =  CHI O/RTODEG  -  PI/2. DO 

c 

T  -  TO 

ALPHAO  =  ALPHAO/RTODEG 
ALPHA  =  ALPHAO 
MUO  =  MUO/RTODEG 
MU  =  MUO 
C 

IF  (WTFUEL  .GT.  WTFLMX )  THEN 
WRITE (6,*)  'WARNING:' 

WRITE( 6 , * )  '  THE  FUEL  WEIGHT  IS  GREATER  THAN  ALLOWED  !!' 

WRITE ( 6 , * )  '  IT  HAS  BEEN  RESET  TO  THE  MAXIMUM  FUEL  WEIGHT.' 

WRITE (9,*)  'WARNING:' 


B-ll 


o  o  o 


c 


c 


c 

c 


c 


c 


WRITE (  9  ,  *  )  '  THE  FUEL  WEIGHT  IS  GREATER  THAN  ALLOWED  !!' 
WRITE ( 9 , * )  '  IT  HAS  BEEN  RESET  TO  THE  MAXIMUM  FUEL  WEIGHT.' 

WTFUEL  =  WTFLMX 
END  IF 


IF  (WTFUEL  .LT.  BINGO)  THEN 


WRITE (  6  ,  *  ) 
WRITE { 6 , * ) 
WRITE (9,*) 
WRITE ( 9 , * ) 
END  IF 


'WARNING:  ' 

'  THE  FUEL  WEIGHT  IS  LESS  THAN 
'WARNING: ' 

'  THE  FUEL  WEIGHT  IS  LESS  THAN 


THE  BINGO  FUEL 
THE  BINGO  FUEL 


i  i  ' 
i  j  ' 


IF  (WTFUEL  .LT 
WRITE (6,*) 
WRITE (6,*) 
WRITE ( 9 , * ) 
WRITE (9,*) 
ENDIF 


.  0.D0)  THEN 
'WARNING: ' 

'  FUEL  WEIGHT  IS  LESS  THAN 
'WARNING: ' 

'  FUEL  WEIGHT  IS  LESS  THAN 


ZERO  ! ! ' 
ZERO  ! 1 ' 


ACM  =  ( WTAC  +  WTFUEL  +  WTARM ) /G 
H  =  -Z 

ROH  =  RHO(H) 

MACH  =  V/DSQRT (1.4D0*P(H) /ROH ) 

CALL  SPLINT ( ALT, SPDMAX, SPMXSD , NMCHMX, H , MXMACH ) 

IF  (MACH  .GT.  MXMACH)  THEN 
WRITE (6,*)  'WARNING:' 

WRITE ( 6 , * )  '  INITIAL  SPEED  EXCEEDS  MAXIMUM  VELOCITY  !!' 

WRITE (6,*)  '  VALUE  RESET  TO  MAXIMUM  VELOCITY.' 

WRITE) 9,*)  'WARNING:' 

WR I T E ( 9 , * )  '  INITIAL  SPEED  EXCEEDS  MAXIMUM  VELOCITY  !)' 
WRITE ( 9 , * )  '  VALUE  RESET  TO  MAXIMUM  VELOCITY.' 

MACH  =  MXMACH 

V  =  MACH  *  DSQRT(1.4DO*P(H)/ROH) 

ENDIF 


MG  =  ACM*G 
VSQ  =  V**2 

REFOR  =  0 . 5DO*ROH*VSQ*AREA 
SNG  -  DSIN(GM) 

CALL  SPLIN2 ( AOACL , MACHCL , CLI FT , CLFTSD , NAOACL , NMCHCL , 

&  ALPHA, MACH, CL) 

CALL  SPLIN2 ( CLFTCD , MACHCD , CDRAG , CDRGSD , NCLCD , NMCHCD , CL , MACH , CD ) 
TH  =  ( MG*SNG  +  REFOR*CD)/DCOS( ALPHA) 

XI  =  T 
X2  =  X 
X3  =  Y 
X4  =  Z 
X5  =  V 
X6  =  GM 
X7  =  CHI 
X8  =  ALPHA 
X9  =  MU 

CALL  ATTITUDE ( Xl , X2 , X3 , X4 , X5 , X6 , X7 , X8 , X9 ) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
COMMANDS  ARE  READ  IN  AND  THE  APPROPRIATE  SUBROUTINES  ARE  CALLED. 
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onnoonnooo 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


FLAG  =  .TRUE. 

DO  WHILE  (FLAG) 

READ( 2,10,ERR=999)  ORDERS 


IF  (ORDERS  .EQ. 

'WA' 

)  THEN 

CALL 

WAYPT ( X , 

V,Z, 

V , GM , CHI , T ) 

ELSE  IF 

(ORDERS 

•  EQ. 

'HE'  ) 

THEN 

CALL 

HEADIN 

ELSE  IF 

( ORDERS 

•  EQ. 

'PI '  ) 

THEN 

CALL 

PITCH 

ELSE  IF 

(ORDERS 

•  EQ. 

'HP'  ) 

THEN 

CALL 

HDGPTC 

ELSE  IF 

(ORDERS 

•  EQ. 

'DR'  ) 

THEN 

CALL 

DRPSTR 

ELSE  IF 

(ORDERS 

.EQ. 

'TF'  ) 

THEN 

CALL 

TERFOL 

ELSE  IF 

(ORDERS 

•  EQ. 

'EN'  ) 

THEN 

FLAG 

=  .FALSE 

• 

ELSE 

WRITE ( 6 / * )  '  COMMAND  ORDERS,'  DOESN '  '  T  EXIST  !!!' 
WRITE ( 9 , * )  '  COMMAND  ' , ORDERS , '  DOESN' 'T  EXIST  !!!' 
FLAG  -  .FALSE. 

END  IF 
END  DO 
C 

999  CONTINUE 
RETURN 

C  *  *  * 

5  FORMAT( 4F10 . 0 ) 

10  FORMAT ( A ) 

C  *** 

END 
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SUBROUTINE  ATTITUDE ( XI ,X2 , X3 , X4 ,X5 , X6 , X7 , X8 , X9 ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C  c 
C  SUBROUTINE  ATTITUDE  IS  CALLED  ONCE  FROM  SUBROUTINE  CONTROL.  IT  IS  C 
C  CALLED  THEREAFTER  AT  THE  END  OF  THE  RK  SUBROUTINE.  ITS  FUNCTION  IS  C 
C  TO  COMPUTE  THE  BODY  ATTITUDE  ANGLES  OF  THE  AIRCRAFT  RELATIVE  TO  THE  C 
C  INERTIAL  REFERENCE  FRAME.  IT  THEN  DUMPS  OUT  THE  AIRCRAFT'S  POSITION  C 
C  SPEED  AND  BODY  ATTITUDE  EVERY  SECOND.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  ' CMNBLK/LIST' 

C 

REFOR  =  . 5DO*RHO ( -X4 ) * ( X5**2 ) *AREA 
C 

ALNG  -  ( TH*DCOS ( X8 ) -REFOR*CD ) /ACM 
ALNG  =  ALNG— G*DSIN(X6) 

C 

FL  =  REFOR*CL 

APTCH  =  ( ( FL+TH*DSIN(X8 ) ) *DCOS < X9 ) -ACM*G*DCOS ( X6 ) )/(ACM*X5) 

C 

AYW  =  ( (FL+TH*DSIN(X8) )*DSIN(X9) )/< ACM*X5*DCOS( X6 ) ) 

C 

ALNG  =  ALNG/G 

APTCH  =  APTCH*X5/G 

AYW  =  AYW*X5*DCOS(X6 )/G 

ATOTAL  =  DSQRT ( ALNG*  *2+APTCH*  *  2+AYW**2 ) 

X2  =  X2/MTOFT 
X3  =  — X3/MTOFT 
X4  =  -X4/MTOFT 
X5  =  X5/MT0FT 
X6  *  X6 
X7  =  X7 
X8  =  X8 

CAL  =  DCOS ( X8  > 

SAL  -  DSIN(X8) 

CGM  =  DCOS ( X6 ) 

SGM  =  DSIN(X6) 

CCH  =  DCOS ( X7 ) 

SCH  =  DSIN(X7) 

CMU  =  DCOS ( X9 ) 

SMU  =  DSIN(X9) 

AA  =  CAL*CGM*CCH-SAL* ( CMU*SGM*CCH+SMU*SCH ) 

BB  =  CAL*CGM*SCH— SAL* ( CMU*SGM*SCH— SMU*CCH ) 

CC  =  ( — CAL*SGM— SAL*CMU*CGM ) 

BCHI  =  DATAN ( BB/AA ) *RTODEG 

IF  ((AA  .LT.  0 . DO )  .AND.  (BB  .LE.  O.DO))  BCHI=180 .DO+BCHI 
IF  ((AA  .LT.  O.DO)  .AND.  (BB  .GT.  O.DO))  BCHI=180 .DO+BCHI 
BCHI  «  90. DO+BCHI 

BGM  =  DATAN ( — CC/DSQRT ( AA* * 2+BB* *2 ) ) 

BMU  =  ( CGM*SMU/DCOS ( BGM ) ) *RTODEG 

BGM  =  BGM*RTODEG 

X6  =  X6*RTODEG 

X7  =  X7*RTODEG 

X7  =  X7+90.D0 

X8  =  X8*RTODEG 

X9  =  X9*RTODEG 

WRITE ( 7 ,3000)Xl,X2,X3,X4,X5, BCHI , BGM , BMU 
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WRITE (8,4000) XI , ACM , X6 , X7 , X8 , X9 , TH , ALNG , APTCH , AYW , ATOTAL 
RETURN 

C  *** 

3000  FORMAT ( IX , 5F1 1 . 2 , 3F11 . 4 ) 

4000  FORMAT ( 1X,11F10.2) 

C  *** 

END 
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SUBROUTINE  WAYPT ( X ,  Y ,  Z  ,  V ,  GM ,  CHI  ,  T ) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  SUBROUTINE  WAYPOINT  IS  CALLED  FROM  CONTROL  TO  CONSTRUCT  A  FLIGHT  C 
C  SEGMENT.  IT  IS  PASSED  THE  PRESENT  POSITION  OF  THE  AIRCRAFT  AND  THE  C 
C  DIRECTION  OF  THE  VELOCITY  VECTOR.  THE  SUBROUTINE  THEN  READS  IN  ITS  C 
C  COMMANDED  POSITION  AND  CALCULATES  THE  DIRECTIONAL  CHANGE  NECESSARY  TO  C 
C  POINT  THE  VELOCITY  VECTOR  TO  THE  DESIRED  WAYPOINT.  FINALLY,  IT  C 
C  CONSTRUCTS  A  FLIGHT  PATH  TO  ACHIEVE  THE  DESIRED  POSITION.  C 
C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


INCLUDE  ' CMNBLK/LI ST' 

C 

REAL *8  MXMACH 
C 

LOGICAL  FLAG 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  READ  IN  THE  DESTINATION  WAYPOINT  AND  CONVERT  TO  MIELE'S  C 
C  COORDINATE  SYSTEM.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

READ( 2,10)  XF , YF , ZF , VF , GEES 
READ ( 2 , * ) 

C 

XF  =  XF*MTOFT 
YF  =  -YF*MTOFT 
ZF  =  -ZF*MTOFT 
VF  «  VF*MTQFT 
C 

H  a  -ZF 
ROH  =  RHO(H) 

MACH  =  VF/DSQRT(1.4D0*P(H)/ROH) 

CALL  SPLINT ( ALT , SPDMAX , SPMXSD , NMCHMX , H , MXMACH ) 

IF  (MACH  .GT.  MXMACH)  THEN 
WRITE ( 6 , * )  'WARNING:' 

WRITE ( 6 , * )  '  COMMANDED  SPEED  EXCEEDS  MAXIMUM  VELOCITY  !!' 

WRITE( 6 , * )  '  VALUE  RESET  TO  MAXIMUM  VELOCITY.' 

WRITE! 9,*)  'WARNING:' 

WRITE ( 9 , * )  '  COMMANDED  SPEED  EXCEEDS  MAXIMUM  VELOCITY  I!' 

WRITE ( 9 , * )  '  VALUE  RESET  TO  MAXIMUM  VELOCITY.' 

MACH  =  MXMACH 

VF  =  MACH  *  DSQRT (1.4D0*P(H) /ROH ) 

END  IF 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  BEGIN  A  DO  LOOP  TO  COMMAND  THE  VELOCITY  VECTOR  INTO  THE  DIRECTION  C 
C  OF  THE  WAYPOINT.  INTEGRATION  IS  PERFORMED  UNTIL  THE  DIFFERENCE  C 
C  BETWEEN  THE  COMMANDED  VELOCITY  VECTOR  AND  THE  CURRENT  VECTOR  IS  C 
C  MINIMAL.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

FLAG  =  .TRUE. 

DO  WHILE  (FLAG) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  COMPUTE  THE  RANGE  AND  THE  GROUND  RANGE  BETWEEN  THE  AIRCRAFT'S  C 
C  PRESENT  POSITION  AND  ITS  DESTINATION.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C 


C 


c 


c 


DELX  =  XF  -  X 
DELY  =  YF  -  Y 
DELZ  =  ZF  -  Z 

VX  =  V*DCOS ( GM ) *DCOS ( CHI ) 

VY  -  V*DCOS ( GM ) *DSIN( CHI ) 

VZ  =*  — V*DS  IN  (  GM ) 

RANGE  =  DSQRT(DELX**2  +  DELY**2  +  DELZ**2) 

OMGALOSX  =  — ( DELY*VZ-DELZ* VY ) /RANGE*  *2 
OMGALOSY  =  -( DELZ *VX-DELX*VZ ) /RANGE* *2 
OMGALOSZ  =■  -  (  DELX*VY-DELY*VX  )  /RANGE**  2 

DELV  =  VF  -  V 
AVEV  =  (VF  +  V)/2.D0 
TTI  =  RANGE/AVEV 
VDOT  =  DELV/TTI 

DRAGACCEL  =  ( CD*REFOR/ACM ) +G*DSIN ( GM ) 

IF  ((VDOT  .LT.  0.D0)  .AND.  ( -VDOT  . GT .  DRAGACCEL))  THEN 
VDOT  =  DSIGN ( DRAGACCEL , VDOT ) 

END  IF 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  THE  PITCH  AND  HEADING  CHANGES  THAT  WILL  BE  NEEDED  TO  ACHIEVE  THE  C 
C  WAYPOINT  ARE  CALCULATED.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


IF  (DSIN(CHI)  .EQ.  0.D0)  THEN 

GMDOT  =  KAY* ( OMGALOSY— OMGALOSZ *DTAN ( GM ) *DSIN ( CHI ) ) /DCOS ( CHI ) 

ELSE 

GMDOT  -  KAY* ( OMGALOSZ *DTAN( GM ) *DCOS ( CHI ) -OMGALOSX ) /DSIN ( CHI ) 

ENDIF 

C 

CHIDOT  =  KAY  *  OMGALOSZ/ (( DCOS ( GM ))* *2 ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  c 

c  COMPUTE  YAW  AND  PITCH  ACCELERATIONS  AT  THESE  COMMANDED  RATES.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

RATIO  =  DSQRT(  ( GMDOT**  2-t- (  CHIDOT*DCOS  (  GM )  )**2)*V**2 
1  +VDOT**2 )/( GEES*G ) 

IF  (RATIO  .GT.  1.0D0)  THEN 
CHIDOT  =  CHIDOT/RATIO 
GMDOT  =  GMDOT/RATIO 
VDOT  =  VDOT/RATIO 
ENDIF 
C 
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non 


CHIDOTI  =  CHIDOT 

CHIDOT  =  F CHI DOT ( CHIDOTI ) 

GMDOTI  =  GMDOT 

GMDOT  =  FGMDOT ( GMDOTI ) 

VDOTI  =  VDOT 

VDOT  =  FVDOT( VDOTI) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTEGRATION  IS  PERFORMED  IF  THE  DIFFERENCES  WARRANT  IT. 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
TF  =  T  +  DELTA 

DO  WHILE  ( (NINT(T*COUNT)  .LT.  NINT( TF*COUNT ) )  .AND.  FLAG) 

CALL  AERO(X, Y,Z,V,GM, CHI, T) 

CALL  RK{ X, Y , Z , V, GM, CHI , T ) 

DELX  =  XF  -  X 
DELY  =  YF  -  Y 
DELZ  =■  ZF  -  Z 

RANGE  =  DSQRT ( DELX*  *2  +  DELY**2  +  DELZ**2) 

IF  (RANGE  .LT.  328. DO)  FLAG= . FALSE. 

END  DO 
END  DO 
RETURN 

C  *  *  * 

10  FORMAT( 5F10 . 0 ) 

C  *  *  * 

END 
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non 


DOUBLE  PRECISION  FUNCTION  FALPHAf ALPHAI ) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

c  FUNCTION  FALPHA  IS  CALLED  FROM  SUBROUTINE  AERO.  THIS  FUNCTION'S  C 
C  JOB  IS  TO  FILTER  THE  ANGLE  OF  ATTACK.  A  1ST— ORDER  LINEAR  FILTER  IS  C 
C  USED  IN  THE  CALCULATION.  ALP HA I ,  THE  PASSED  PARAMETER,  IS  THE  C 
C  COMMANDED  ANGLE  OF  ATTACK.  ALPHAO  IS  THE  PRESENT  HELD  VALUE  OF  ANGLE  C 
C  OF  ATTACK.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  'CMNBLK/LIST' 

C 

TERM  =  DEXP ( -DT/AOATAU ) 

FALPHA  =  ALPHAO *TERM  +  ALPHAI *( 1 . DO-TERM ) 

ALPHAO  =  FALPHA 

RETURN 

END 

CCCCC 

CCCCC 

CCCCC 

DOUBLE  PRECISION  FUNCTION  FMU(MUI) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FUNCTION  FMU  IS  CALLED  FROM  SUBROUTINE  AERO.  THIS  FUNCTIONS'S  C 
C  JOB  IS  TO  FILTER  THE  ROLL  ANGLE.  A  1ST-ORDER  LINEAR  FILTER  IS  USED  C 
C  IN  THE  CALCULATION.  MUI ,  THE  PASSED  PARAMETER,  IS  THE  COMMANDED  ROLL  C 
C  ANGLE.  MUO  IS  THE  PRESENT  HELD  ROLL  ANGLE.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  'CMNBLK/LIST' 

C 

REAL *8  MUI 
C 

TERM  -  DEXP ( — DT/MUTAU ) 

FMU  =  MU0*TERM  +  MUI *( 1 . DO-TERM ) 

MUO  =  FMU 

RETURN 

END 

CCCCC 

CCCCC 

CCCCC 

DOUBLE  PRECISION  FUNCTION  FGMDOT ( GMDOTI ) 

CC CCCCCCC CCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCC CCCCC CCCCC 


c  c 

C  FUNCTION  FGMDOT  IS  CALLED  FROM  VARIOUS  SUBROUTINES.  THIS  C 

C  FUNCTIONS'S  JOB  IS  TO  FILTER  THE  PITCH  ACCELERATION.  A  1ST-ORDER  C 

C  LINEAR  FILTER  IS  USED  IN  THE  CALCULATION.  GMDOTI,  THE  PASSED  C 

C  PARAMETER,  IS  THE  COMMANDED  PITCH  ANGLE.  GMDOTO  IS  THE  PRESENT  HELD  C 

c  pitch  Angle.  c 

c  c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INCLUDE  'CMNBLK/LIST' 

c 

TERM  =  DEXP ( -DELTA/GMTAU } 
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FGMDOT  =  GMDOTO  *TERM  +  GMDOTI *( 1 . DO-TERM ) 

GMDOTO  =  FGMDOT 

RETURN 

END 

CCCCC 

CCCCC 

CCCCC 

DOUBLE  PRECISION  FUNCTION  FCHIDOT( CHIDOTI ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FUNCTION  FCHIDOT  IS  CALLED  FROM  VARIOUS  SUBROUTINES.  THIS  C 
C  FUNCTIONS'S  JOB  IS  TO  FILTER  THE  YAW  ACCELERATION.  A  1ST-ORDER  C 
C  LINEAR  FILTER  IS  USED  IN  THE  CALCULATION.  CHIDOTI,  THE  PASSED  C 
C  PARAMETER,  IS  THE  COMMANDED  YAW  ANGLE.  CHIDOTO  IS  THE  PRESENT  C 
C  HELD  YAW  ANGLE.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
INCLUDE  ' CMNBLK/LIST ' 

c 

TERM  =  DEXP ( — DELTA/CHITAU ) 

FCHIDOT  =  CHIDOTO* TERM  +  CHIDOTI *( 1 . DO-TERM ) 

CHIDOTO  -  FCHIDOT 

RETURN 

END 

CCCCC 

CCCCC 

CCCCC 

DOUBLE  PRECISION  FUNCTION  FVDOT(VDOTI) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FUNCTION  FVDOT  IS  CALLED  FROM  VARIOUS  SUBROUTINES.  THIS  C 
C  FUNCTIONS 'S  JOB  IS  TO  FILTER  THE  LONGITUDNAL  ACCELERATION.  A  C 
C  1ST-ORDER  LINEAR  FILTER  IS  USED  IN  THE  CALCULATION.  VDOTI , THE  C 
C  PASSED  PARAMETER,  IS  THE  COMMANDED  ACCELERATION.  VDOTO  IS  THE  PRESENT  C 
C  HELD  ACCELERATION.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
INCLUDE  'CMNBLK/LIST' 

c 

TERM  =  DEXP( -DELTA/VTAU ) 

FVDOT  =  VDOTO *TERM  +  VDOTI *( 1 . DO-TERM ) 

VDOTO  =  FVDOT 

RETURN 

END 

CCCCC 

CCCCC 

CCCCC 

DOUBLE  PRECISION  FUNCTION  RHO(H) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  FUNCTION  RHO  COMPUTES  THE  STANDARD  ATMOSPHERIC  DENSITY  AS  A  C 
C  FUNCTION  OF  ELEVATION,  H.  THE  PASSED.  PARAMETER  H  HAS  UNITS  FT.  C 
C  THE  DENSITY  HAS  UNITS  SLUGS/FT**3.  C 
C  C 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


INCLUDE  ' CMNBLK/LIST' 

c 

IF  (H  .GT.  3.6089D4)  THEN 

RHO  =  RHO0*0 . 297lD0*DEXP( ( 3 . 6089D4-H )/2 . 08067D4 ) 

ELSE 

THETA  =  l.D0-(6.875D-6}*H 

RHO  =  RHO0  *DEXP (4.2561DO*DLOG{ THETA )  ) 

END  IF 

RETURN 

END 

CCCCC 

CCCCC 

CCCCC 

DOUBLE  PRECISION  FUNCTION  P(H) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  FUNCTION  P ( H )  COMPUTES  THE  STANDARD  PRESSURE  AS  A  FUNCTION  OF 

C  ELEVATION.  THE  PASSED  PARAMETER  H  HAS  UNITS  FT.  THE  PRESSURE 

C  HAS  UNITS  LBS/FT* *2 . 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  'CMNBLK/LIST' 

c 

IF  (H  .GT.  3.6089D4)  THEN 

P  =  P0*0 . 223  4D0*DEXP ( ( 3 . 6  089D4-H ) /2 . 08  067D4 ) 

ELSE 

THETA  =  l.D0-(6.875D-6)*H 
P  =  P0*DEXP( 5 . 256lDO*DLOG( THETA) ) 

END  IF 

RETURN 

END 
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SUBROUTINE  AERO ( X , Y , Z , V , GM , CHI , T ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  SUBROUTINE  AERO  IS  CALLED  BY  THE  VARIOUS  CONTROL  ROUTINES.  ITS  C 
C  FUNCTION  IS  TO  CALCULATE  THE  AERODYNAMIC  FORCES  BASED  ON  THE  C 
C  COMMANDED  ACCELERATIONS.  THE  CONTROL  VARIABLES  CALCULATED  ARE  THE  C 
C  ANGLE  OF  ATTACK  AND  THE  VELOCITY  ROLL  ANGLE.  THE  THRUST  NEEDED  IS  ALSO  C 
C  CALCULATED.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC.CCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  ' CMNBLK/LIST ' 

C 

REAL *8  MU I 
C 

H  =  -Z 

ROH  =>  RHO(H) 

MACH  =  V/DSQRT ( 1 . 4  DO  *  P ( H ) /ROH ) 

MG  =  ACM*G 
MV  =  ACM*V 
VSQ  =  V**2 

REFOR  =  0 . 5DO*ROH*VSQ*AREA 
SNG  =  DSIN(GM) 

CSG  =  DCOS(GM) 

SNCHI  =  DSIN(CHI) 

CSCHI  =  DCOS(CHI) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  COMPUTE  THE  ROLL  ANGLE  FROM  THE  COMMANDED  PITCH  AND  YAW  RATES.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

MU I  =  ATAN2{ (CHIDOT*CSG) , (GMDOT+G*CSG/V) ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  FILTER  THE  ROLL  ANGLE.  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

MU  =  FMU(MUI) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CALL  NEWTON  TO  DETERMINE  THE  ANGLE  OF  ATTACK  CORRESPONDING  TO  THE  C 
C  FILTERED  ROLL  ANGLE,  COMMANDED  PITCH  RATE  AND  LONGITUDINAL  ACCELERATION . C 
C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


ALPHAI  =  ALPHA 
CALL  NWTN1D( ALPHAI ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  INTERPOLATE  TO  FIND  THE  MAXIMUM  ANGLE  OF  ATTACK  AT  THIS  MACH  C 
C  NUMBER.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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IF  ( DABS ( ALPHAI )  .GT.  2.D0*PI)  THEN 
ALPHAI  =  DMOD { ALPHA I , 2 . DO  *  PI ) 

END  IF 

IF  ( DABS ( ALPHAI )  .GT.  PI)  THEN 

TEMALF  =  2 . D0*PI  -  DABS (ALPHAI ) 

ALPHAI  =  DSIGN( TEMALF, -ALPHAI ) 

END  IF 

CALL  SPLINT ( AIRSPD , AOAMAX , AAMXSD , NAOAMX , MACH , AL FMAX ) 
IF  (ALPHAI  .GT.  ALFMAX)  THEN 
ALPHAI  =  ALFMAX 
ENDIF 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  FILTER  THE  ANGLE  OF  ATTACK.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

ALPHA  =  FALPHA( ALPHAI ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  DETERMINE  THE  LIFT  COEFFICIENT  AT  THIS  VALUE  OF  ANGLE  OF  ATTACK.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  SPLIN2 ( AOACL , MACHCL , CL I  FT , CLFTSD , NAOACL , NMCHCL , 

&  ALPHA, MACH, CL) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  DETERMINE  THE  DRAG  COEFFICIENT  AT  THIS  LIFT  COEFFICIENT.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  SPLIN2 ( CLFTCD ,MACHCD , CDRAG , CDRGSD , NCLCD , NMCHCD , CL , MACH , CD ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  DETERMINE  THE  MAXIMUM  THRUST  AT  THIS  ALTITUDE  AND  MACH  NUMBER.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  SPLIN2 ( ALTMIL , MCHMIL , PWRMI L , PRMLSD , NALTML , NMCHML , 

&  H , MACH , THMAX ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  COMPUTE  THE  THRUST  FOR  THIS  COMMANDED  LONGITUDINAL  ACCELERATION.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

TH  =  ( ACM*VDOT+MG*SNG+CD*REFOR )/DCOS ( ALPHA ) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 
c  LIMIT  THE  THRUST.  C 
C  C 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


IF  (TH  .GT.  THMAX)  THEN 
TH=THMAX 
END  IF 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMPUTE  THE  PERCENT  MAXIMUM  RPMs ,  AND  DETERMINE  THE  THRUST 
SPECIFIC  FUEL  CONSUMPTION. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

RPM  =  ( TH/THMAX) **ROOT 

CALL  SPLINT ( PCTRPM , TSFC , TS FCSD , NCTRPM , RPM , CT ) 

C 

RETURN 

END 
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noon 


SUBROUTINE  NWTNID(X) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  NEWTON  PERFORMS  A  1-DIMENSIONAL  ROOT  FINDING  ALGORITHM 
USING  THE  NEWTON-RAPHSON  METHOD. 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IMPLICIT  REAL *8  (A-H,0-Z) 

C 

LOGICAL  FLAG 
C 

EPSLNX  =  . 00875D0 
FLAG  =  .TRUE. 

JCOUNT  =  0 
DO  WHILE  (FLAG) 

C 

FUNC1  =  FUNCFl(X) 

C 

DF1X  »  FDFlX(X) 

C 

DET  =  DF1X 

C 

DETINV  =  1.D0/DET 

C 

XNEW  -  X  -  DETINV*FUNC1 

C 

IF  ( DABS ( X— XNEW )  .LE.  EPSLNX)  THEN 
FLAG  =  .FALSE. 

ENDIF 

C 

JCOUNT  -  JCOUNT  +  1 
IF  (JCOUNT  .EQ.  50)  THEN 
XNEW  =>  (XNEW  +  X)  /2  .  DO 
ELSE  IF  (JCOUNT  .EQ.  1000)  THEN 

WRITE (6,*)  '  NEWTON  OVERLOAD  l!' 

WRITE (9,*)  '  NEWTON  OVERLOAD  II' 

STOP  '  MAJOR  ERROR  ! ! ' 

ENDIF 

C 

X  =  XNEW 
C 

END  DO 
C 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  FUNCFl ( ALPHAI ) 

C 

INCLUDE  'CMNBLK/LIST' 

C 

CALL  SPLIN2 ( AO ACL , MACHCL , CLI FT , CLFTSD , NAOACL , NMCHCL , 

&  ALPHAI , MACH, CL) 

CALL  SPLIN2 ( CLFTCD , MACHCD , CDRAG, CDRGSD , NCLCD , NMCHCD , CL , MACH , CD ) 
C 

FUNCFl  =  ((( (ACM*VDOT+MG*SNG+REFOR*CD )*DTAN (ALPHAI )+REFOR*CL) 

&  *DCOS ( MU )  -  MG*CSG)  /MV)  -  GMDOT 

RETURN 
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END 


C 

C 


C 


DOUBLE  PRECISION  FUNCTION  FDFlX( ALPHA I ) 
INCLUDE  ' CMNBLK/LIST ' 


CALL  SPLIN2 ( AOACL , MACHCL , DCLDA , DCLDASD , NAOACL , NMCHCL , 

&  ALPHAI , MACH , DERCL ) 

CALL  SPLIN2 ( CLFTCD , MACHCD , DCDDCL , DCDDCLSD , NCLCD , NMCHCD , 
&  CL , MACH , DERCD ) 


FDFlX  = 

& 

& 

FDFlX  = 

& 

RETURN 

END 


(  RE FOR* CD* ( 1 . DO/DCOS ( ALPHAI )  )  ** 2  + 
REFOR*DERCD*DERCL*DTAN( ALPHAI )  + 

REFOR* DERCL  )  *  DCOS(MU)  /  MV 
( ( ACM*VDOT  +  MG*SNG)*( 1. DO/DCOS (ALPHAI ) )**2) 
*(DCOS(MU)/MV)  +  FDFlX 
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SUBROUTINE  RK ( X, Y , Z , V, GM , CHI , T ) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  RK  IS  CALLED  FROM  THE  VARIOUS  CONTROL  ROUTINES.  ITS 
FUNCTION  IS  TO  IMPLEMENT  THE  4TH-ORDER  RUNGE-KUTTA  FORMULA.  SUBROUTINE 
RK  ADVANCES  A  TIME-STEPPED  SOLUTION  TO  THE  DIFFERENTIAL  SYSTEM 
DESCRIBING  THE  TIME  RATE  OF  CHANGE  OF  THE  AIRCRAFT'S  INERTIAL  POSITION 
AND  VELOCITY,  AND  MASS. 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c 

c 


c 


c 


INCLUDE  ' CMNBLK/LIST ' 

CALL  F ( X , Y , Z , V , GM , CHI , ACM , T ) 

Fll  =  Fl 
F21  =  F2 
F31  =  F3 
F41  »  F4 
F51  =  F5 
F61  =  F6 
F7 1  =  F7 

XSIN  =  X+DT*0 . 5D0*Fll 
YSIN  =  Y+DT*0.5DO*F21 
ZSIN  =  Z+DT*0 . 5D0*F31 
VSIN  =  V+DT*0 . 5D0*F41 
GMSIN  =  GM+DT*0.5D0*F51 
CHSIN  -  CHI+DT*  0 . 5D0*F61 
ACMSIN  =*  ACM+DT*0.5D0*F71 
T  =  T+DT*  0 . 5D0 

CALL  F ( XSIN, YSIN, ZSIN, VSIN , GMSIN, CHSIN, ACMSIN, T) 

F12  =  FI 
F22  =  F2 
F32  *=  F3 
F42  =■  F4 
F52  -  F5 
F62  -  F6 
F72  -  F7 

XSSIN  =  X+DT*0 . 5D0*F12 
YSSIN  -  Y+DT*0 . 5D0*F22 
ZSSIN  =  Z+DT*0 . 5D0*F32 
VSSIN  -  V+DT*0.5D0*F42 
GMSSIN  =  GM+DT*0.5D0*F52 
CHSSIN  =  CHI+DT*0 . 5D0*F62 
ACMSSIN  =  ACM+DT*  0 . 5DQ  *  F7  2 

CALL  F( XSSIN, YSSIN, ZSSIN, VSSIN, GMSSIN, CHSSIN, ACMSSIN, T) 


F13  = 

Fl 

F23  = 

F2 

F33  - 

F3 

F43  - 

F4 

F5  3  = 

F5 

F63  = 

F6 

F73  = 

F7 

XSPl 

=  X+DT*Fl3 

YSP 1 

=  Y-t-DT*F2  3 

ZSPI 

=  Z+DT*F33 

VSP1 

=  V+DT*F43 
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GMSPl  =  GM+DT*F53 
CHSP1  =  CHI +DT*F6  3 
ACMSP1  =  ACM+DT*F73 
T  =  T+DT*  0 . 5D0 

CALL  F ( XSPl ,YSP1,ZSP1,VSP1, GMSPl , CHSP1 , ACMSPl , T ) 

F14  =  Fl 
F24  =  F2 
F34  =  F3 
F44  =  F4 
F5  4  =  F5 
F6  4  =  F6 
F74  =  F7 

X  =  X+DT* (F11/6.D0+F12/3. D0+F1 3/3 . DO+Fl 4/6 . DO ) 

Y  =  Y+DT* (F21/6.DO+F22/3.DO+F23/3.DO+F24/6.DO) 

Z  =  Z+DT* (F31/6.DO+F32/3.DO+F33/3.DO+F34/6.DO) 

V  =  V+DT* ( F41/6 .D0  +  F42/3.D0  +  F43/3. D0+F4  4/6 . DO ) 

GM  =  GM+DT* ( F51/6 . DO+F52/3 . D0+F53/3 .DO+F54/6.DO) 
CHI  =  CHI+DT* (F61/6.D0+F62/3.D0+F63/3.D0+F64/6.D0) 
ACM  =  ACM+DT* (F71/6.O0+F72/3.D0+F73/3.D0+F74/6.D0) 
IF(CHI.LT.-PI)  CHI=CHI+2.D0*PI 

IF(CHI.GT.PI)  CHI=CHI-2 .DO* PI 

JFLAG  =  NINT(T*COUNT) 

IF  (MOD( JFLAG, ICOUNT)  .EQ.  0)  THEN 
XI  =  T 
X2  =  X 
X3  =  Y 
X4  =  Z 
X5  =  V 
X6  =  GM 
X7  =  CHI 
X8  =  ALPHA 
X9  =  MU 

CALL  ATTITUDE ( XI , X2 , X3 , X4 , X5 , X6 , X7 , X8 , X9 ) 

ENDIF 

RETURN 

END 
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SUBROUTINE  F ( X , Y , Z , V , GM , CHI , MAC , T ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  SUBROUTINE  F  COMPUTES  THE  INSTANTANEOUS  TIME  RATES  OF  CHANGE  OF  C 
C  THE  VEHICLE'S  CARTESIAN  POSITION,  SPEED,  VELOCITY  PITCH  AND  YAW  ANGLES  C 
C  AND  MASS.  SUBROUTINE  F  IS  CALLED  BY  SUBROUTINE  RK .  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INCLUDE  ' CMNBLK/LIST' 

C 

REAL *8  MAC 
C 

H  =  -Z 

ROH  *  RHO(H) 

MG  =  MAC*G 
MV  =  MAC*V 
VSQ  »  V**2 

REFOR  =  0 . 5DO*ROH*VSQ*AREA 
SNG  -  DSIN(GM) 

CSG  -  DCOS(GM) 

SNCHI  =  DS IN ( CHI ) 

CSCHI  =  DCOS(CHI) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


c  c 

C  COMPUTE  THE  COMPONENT  OF  VELOCITY  ALONG  THE  X-AXIS  IN  MIELE'S  C 
C  FLAT  EARTH  COORDINATE  SYSTEM.  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
FI  =  V*CSG*CSCHI 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  c 
c  COMPUTE  THE-  COMPONENT  OF  VELOCITY  ALONG  THE  Y-AXIS  IN  MIELE'S  C 
C  FLAT  EARTH  COORDINATE  SYSTEM.  C 
C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


F2 


V*CSG*SNCHI 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  COMPUTE  THE  COMPONENT  OF  VELOCITY  ALONG  THE  Z-AXIS  IN  MIELE'S  C 

C  FLAT  EARTH  COORDINATE  SYSTEM.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

F3  =  -V*SNG 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  '  c 

c  COMPUTE  THE  LONGITUDINAL  ACCELERATION  V-DOT.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


F4  =  (TH*DCOS( ALPHA )-REFOR*CD) /MAC 
F4  =  F4-G*SNG 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  COMPUTE  THE  RATE  OF  CHANGE  OF  THE  VELOCITY  PITCH  ANGLE  GAMMA-DOT.  C 
C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


FL  =  REFOR*CL 

F5  =  ( ( FL+TH*DS IN( ALPHA ) ) *DCOS { MU ) -MAC*G*CSG ) /( MAC* V ) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  COMPUTE  THE  RATE  OF  CHANGE  OF  THE  VELOCITY  YAW  ANGLE  CHI-DOT.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

F6  =  ( ( FL+TH*DSIN( ALPHA ) ) *DSIN( MU ) )/( MAC*V*CSG ) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  COMPUTE  THE  RATE  OF  CHANGE  OF  THE  MASS.  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


F7  =  -CT*TH/(G*3600.D0) 

c 

ALNG  =  F  4/G 
APTCH  =  F5*V/G 
AYW  =  F6*V*CSG/G 

ATOTAL  =  DSQRT ( ALNG*  *2+APTCH*  *  2+AYW**  2 ) 
IF  (ATOTAL  .GT.  GEES)  THEN 

WRITE ( 6 ,1111 )  T, ALNG, APTCH, AYW, ATOTAL 
WRITE (9,1111)  T, ALNG, APTCH, AYW, ATOTAL 
ENDIF 
RETURN 

c  *  *  * 

1004  FORMAT ( 1H  , ' MIL ' ) 

1005  FORMAT ( 1H  , ' AB ' ) 

1013  FORMAT( IX, 'BINGO  FUEL') 

1111  FORMAT(lH  , ' FT I ME , AL ,AP , AY , AT= ' , 5F10 . 4 ) 

C  *** 

END 
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APPENDIX  C 


MODMOFO.COM 


This  appendix  contains  the  source  code  to  the  file  MODMOFO.COM.  It  also 
contains  a  ".COM"  file  that  establishes  the  common  block  cross-reference  files.  As  alluded 
to  earlier,  the  use  of  these  files  will  make  modifications  to  the  source  code  much  easier. 

MODMOFO.COM  is  a  derivative  of  the  source  code  used  in  modifying  another 
model  within  the  Advanced  Systems  Division.  The  original  code,  MODALARM.COM, 
was  authored  by  Ciro  L.  Pinto-Coelho.  The  user  will  find  the  source  code  internally 
documented  for  immediate  implementation.  The  insightful  user  will  find  it  applicable  to 
other  models  as  well. 
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$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

s 

$ 

$ 

$ 

s 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 


:  This  command  file  mimics  the  UNIX  make  utility  in  modifying 

!  a  particular  group  of  source  and  common  block  files.  This  command 
!  procedure  relies  on  a  specific  directory  structure.  The  top  level 
!  directory  has  the  same  name  as  the  model.  Below  it  are  a  number  of 
!  other  subdirectories  as  follows: 

!  [MODEL] 

J  |  -  Location  of  executable. 

!  [MODEL. CODE] 

i  |  -  Location  of  source  code  and  this  command 

1  j  procedure. 

!  [ MODEL . CODE . CMNBLKS ] 

!  |  -  Location  of  common  block  files. 

!  [ MODEL . CODE . CMNBLKS . DAT ] 

!  |  Location  of  common  block  cross 

!  j  reference  files. 

!  The  procedure  relies  on  the  coding  practice  of  include  statements 

1  to  incorporate  common  blocks  into  a  program.  The  common  blocks 
!  are  defined  by  files  -  whose  prefix  is  the  name  of  the  common  block  - 
!  in  the  subdirectory  [ MODEL . CODE . CMNBLKS ] .  These  files  have  extension 
!  " . CMN " .  Further,  each  common  block  has  a  cross  reference  file  - 

!  located  in  the  subdirectory  [MODEL. CODE . CMNBLKS . DAT ]  -  that  contains 

1  a  list  of  those  source  code  files  that  use  it.  These  files  can  be 
1  generated  with  the  command  procedure  MKDAT.COM  in  the  subdirectory 
!  [MODEL. CODE. CMNBLKS . DAT ] . 

!  This  proceudre  relies  on  the  fact  that  modified  source  or  common 

!  block  files  have  version  numbers  greater  than  one.  Therefore, 

1  this  command  procedure  generates  a  list  of  the  source  code  files  that 
1  have  been  modified  and  places  them  in  a  temporary  file.  It  then 
!  generates  a  list  of  those  common  block  files  that  'have  been  modified. 

!  These  names  are  saved  and  the  cross  reference  file  for  each  common  block 
!  is  read  to  determine  which  source  ocde  files  must  be  compiled.  From  these 
!  two  steps  the  procedure  builds  an  alphabetically  sorted  list  of  the 
!  source  code  files  to  be  recompiled.  The  list  is  sorted  to  avoid 
!  recompiling  a  source  file  more  than  once.  Each  source  file  is  then 
!  recompiled  and  a  copy  is  placed  in  the  text  library  MODEL. TLB,  and  a  copy 

!  of  the  object  file  is  placed  in  the  library  MODEL. OLB.  The  linker  is 

1  then  called  to  create  an  executable  image  name  MODEL.EXE  in  the 
[•subdirectory  [MODEL].  After  each  successful  compilation  the  procedure 
I  purges  the  old  version(s)  of  the  source  code  file  and  revises  its  version 
!  number  to  one.  Recall  that  this  is  the  key  to  the  operation  of  this 

!  procedure.  Finally,  a  bit  of  clean-up  of  the  temporary  file  is  done  before 

i  completion. 


Comments .  .  . 

—  An  option  was  added  to  allow  generation  of  listing  files  from  the 
compilation  step.  The  option  is  activated  by  supplying  an  ”1"  or 
"L",  for  list,  on  the  command-line. 

Establish  error  handling  procedures... 
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S  i 

$  ON  CONTROL  Y  THEN  GOTO  FINISHED 
$  ON  ERROR  THEN  GOTO  FINISHED 
S  ! 

$  I  Define  the  model  name  and  default  file  extension... 

S  ! 

$  CMNDATPTH  =  " D I S K4 : [ MO FO . CODE . CMNBLKS . D AT ] " 

$  CODEDIR  =  "DISK4 : (MOFO. CODE ] " 

$  CMNPTH  =  "DISK4 : [MOFO . CODE . CMNBLKS 1 " 

$  MODEL  =  "MOFO" 

$  EXT  =  "MFO" 

$  ! 

$  IF  PI  .NES.  "  ”  THEN  PI  =  f $ e di t ( P 1 , " LOWE RC AS E * ) 

$  ! 

$  !  Initialize  the  counter  for  files  compiled... 

$  ! 

$  NFILES  =  0 

$  ! 

$  1  Set  the  default  to  the  code  directory,  clear  the  screen  and  begin 
$  1  execution  of  the  procedure... 

$  ! 

$  SET  DEFAULT  'CODEDIR' 

$  SET  TERMI NAL/WI DTH—  8  0 
$  WRITE  SYS  $  OUT  PUT  "  " 

$  WRITE  SYS $OUTP UT  MODEL,"  make  utility..." 

$  WRITE  SYS  $  OUT  PUT  "  " 

$  ! 

$  :  Assign  all-common  block  files  to  logical  names.  This  is  a  model  specific 

$  l  section  of  the  command  procedure. . . 

$  ! 

$  ASSIGN  ' CMNPTH ' CMNBLK . CMN  CMNBLK 

$  ! 

$  !  Generate  a  directory  of  those  model  routines  that  were  modified. 

$  !  Include  those  files  that  are  affected  by  modifications  to 

$  1  their  common  block (s !  .  .  . 

$  ! 

$  WRITE  SYS  $  OUTPUT  "  Creating  directory  of  modified  source  code  files.. 

$  DIRECTORY/VE RSI ON- 1/NOHE ADING/NOTRAILING/OUTPUT® ' MODEL ' . CODE/EXCLUDE* ( * . * 

$  t 

$  WRITE  SYS  $  OUTPUT  "  Creating  directory  of  modified  common  block  files. 

$  DIRE CTORY/ VERS  I ON  =  l/NOHE ADI NG/NOTRAI LI NG/OUTPUT= 'MODEL '  . CMNS/EXCLUDE= (  *  .  * 

$  ! 

$  !  Generate  a  file  that  contains  a  list  of  the  files  affected 

$  !  by  the  common  modification... 

$  ! 

$  WRITE  SYS  $  OUTPUT  "  Creating  composite  directory  listing..." 

$  ! 

$  OPEN/READ  CMNS_FILE  'MODEL'. CMNS 
$  OPEN/WRITE  PURGE_CMNS_FILE  ' MODEL '. PURGE 

$  ! 

$  STARTCMNS : 

$  ! 

$  READ/END _ OF_F I LE=ENDCMNS  CMNS_FILE  MODULE_NAME 

$  MODULE_NAME  =  F $ E DIT ( MODULE_NAME , " COMPRE S S , TRIM " ) 

$  FILE  =  F  $  PARS  E ( MODULE_NAME ,  ,  ,  "NAME "  ) 

$  DFILE  =  CMNDATPTH  +  FILE  +  ".DAT" 


. 1 )  * . 'EXT ' 

n 

. 1 )  [ . CMNBLKS ] * . CMN 
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!  Form  one  file  that  contains  all  the  file  names... 

COPY  *  MODEL  '  . CODE  +  ' DF I LE '  ' MODEL '  . CODE 

i 

CFILE  =  CMNPTH  +  FILE  +  " . CMN " 

LFILE  =  CMNPTH  +  MODEL  +  MCMNn 

t 

WRITE/SYMBOL  PURG£_CMNS  FILE  CFILE 
LIBRARY/REPLACE/TEXT  'LFILE'  'CFILE' 

! 

GOTO  STARTCMNS 
ENDCMNS : 

t 

CLOSE  PURGE_CMNS_FILE 

i  Sort  the  modified  file  names  so  that  the  LAST^FILE  construct 
i  in  the  compilation  loop  will  function  properly... 

\ 

WRITE  SYS$OUTPUT  "  Sorting  composite  directory...” 

SORT  'MODEL'. CODE  'MODEL'. CODE 
PURGE  'MODEL'. CODE 

RENAME  ' MODEL CODE . *  ' MODEL '. CODE . 1 

i 

!  Open  the  directory  file  for  input... 

i 

OPEN/READ  INPUT_FILE  'MODEL'. CODE 

i 

LAST_FILE  =  "  " 

I  Loop  through  the  list  of  modified  routines  and  compile  each,  being 
I  careful  not  to  compile  a  particular  routine  twice... 

i 

BEGIN  : 

i 

READ/END_OF_F I LE=END  INPUT_FILE  MODULE_NAME 
MODULE_NAME  =  F$  EDIT { MODULE_NAME ,  "COMPRESS , TRIM"  ) 

FILE  =  F$PARSE (MODULE_NAME ,7r "NAME" ) 

• 

IF  FILE  .  EQS  .  LAST_FILE  THEN  GOTO  BEGIN 

SOURCE  =  FILE  +  +  EXT 

OBJECT  =  FILE  +  ".OBJ" 

i 

SET  TERMI NAL/WI DTH  =  8  0 
WRITE  SYS  $OUTPUT  "  " 

WRITE  SYS$OUTPUT  "Recompiling  module  -  ", f$edit ( source , "LOWERCASE" ) 
!  Determine  if  a  list  file  is  to  be  created... 

IF  PI  . EQS  .  "1"  THEN  GOTO  LIST 

FORT RAN/ CO NT I NUAT ION=  5  0 /G_F  LOAT ING/CHE  CK«ALL  ' FILE '  .  ' EXT ' 

GOTO  CLEAN 
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LIST : 

FORTRAN/CONTINUATION=50/G_FLOATING/LIST/CHECK*ALL  7  FILE '  .  ' EXT 7 

I 

CLEAN : 

> 

PURGE  7  FILE  7  .  ' EXT  7 

RENAME  7  FILE  7  .  7  EXT  7  .  *  7 F I LE 7  .  7  EXT  7  .  1 

NFILES  «  NF I LE  S  +  1 

WRITE  SYS  $  OUTPUT  "Replacing  text  module  -  " , f$edit (source, " LOWER CAS  E  "  ) 
LIB/REPLACE/TEXT  7  MODEL  '  'FILE'.'EXT7 

i 

WRITE  SYS  $  OUT  P  UT  "Replacing  object  module  -  "  , f  $ e di t  <  object, " LOWERCASE " ) 
LIB/REPLACE/OBJ  7  MODEL  7  'FILE'. OBJ 

DELETE  7  FILE  7  . OBJ . * 

LAS  T^F I LE  =  FILE 

» 

GOTO  BEGIN 

» 

END  : 

t 

J  Compress  the  object  library  to  minimize  the  amount  of  disk 
1  space  used  -  -  . 

j 

LIB/COMPRESS/OB J  'MODEL7 
PURGE  7  MODEL  7  . 0  LB 

RENAME  'MODEL'. OLB  7  MODEL '  . 0 LB . 1 

t 

i  Link  the  model  to  create  an  executable  image  if  compilation  has  proceeded 
l  thi s  point.  .  . 

» 

IF  NFILES  .EQ.  0  THEN  GOTO  FINISHED 

t 

SET  T ERMI NAL/WI DTH  =  8  0 
WRITE  SYS  $  OUTPUT  "  " 

WRITE  S  Y  S  $  OUTPUT  "Relinking  the  " , MODEL , w  model..." 

LINK  7  MODEL  7 /LIB/IN  CLUDE= (  'MODEL7  ) /EXE  CUTABLE  = 'MODEL'  .EXE 

T 

J  Place  the  executable  image  in  its  default  location... 

j 

DELETE  [ - ]  7  MODEL '  . EXE .  * 

RENAME  'MODEL'. EXE  { - ]  7  MO DE L 7  . EXE . 1 


WRITE 

S  YS  $  OUTPUT  "" 

WRITE 

S  YS  $  OUTPUT  MODEL,"  model  modification 

completed  1  n 

WRITE 

S  YS  $  OUTPUT  "" 

Clean- 

up  before  exiting... 

WRITE 

SYS$ OUTPUT  "Purging  common  block  files 

*» 

WRITE 

SYS  $  OUTPUT  "  " 

OPEN/READ  PURGE  CMNS_FILE  7  MODEL '. PURGE 
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BEGINP : 


$ 

$  ! 

S  READ/END  OF  FILE=ENDP  PURGE  CMNS  FILE  MODULE  NAME 

$  ! 

$  PURGE  'MODULE_NAME ' 

$  RENAME  '  MODULE _ NAME '  . *  ' MODULE _ NAME '  . I 

$  GOTO  BEGINP 

$  ! 

$  ENDP : 

$  ! 

$  ! 

$  FINISHED: 

$  ! 

$  !  Close  and  delete  the  directory  file... 

$  ! 

$  CLOSE  INPUT_FILE 
$  CLOSE  CMNS_FILE 
$  CLOSE  PURGE_CMNS_FILE 
$  ! 

$  DELETE  'MODEL' . CODE . * 

$  DELETE  'MODEL' . CMNS . * 

$  DELETE  'MODEL' . PURGE . * 

$  ! 

$  l  Deassign  all  model  specific  logical  symbols... 

$  ! 

$  DEASSIGN  CMNBLK 

$  ! 

$  WRITE  S  Y  S  $  OUTPUT  "Donel" 

$  EXIT 
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l  Command  procedure  to  establish  the  common  block  cross  reference  files... 

ON  CONTROL_Y  THEN  GOTO  FINISHED 
ON  ERROR  THEN  GOTO  FINSHED 


I  Define  some  logicals... 


CMNDAT  =  B  DISK 4  :  [ MOFO . CODE . CMNBLKS . DAT ]  " 
CMNPTH  =  "DISK4 : [MOFO . CODE . CMNBLKS ] » 

CODE  =  "DISK4 : [MOFO .CODE ] * 

MODEL  m  " MOFO " 

CEXT  =  " CMN n 

DEXT  =  "DAT" 

SEXT  =  " MFO " 


I  Generate  a  directory  of  those  model  common  block  files  that  were  modified. 

COMMONS  =  CMNPTH  +  " * . "  +  CEXT 
SFILE  as  CODE  +  +  SEXT 


D I RE CTOR Y/VERS I ON= 1/NOHEAD I NG/NOTRAI LING/OUT PUT** ' MODEL '  . DCT/EXCLUDE= (  *  . *  . 1  > 
» 

I  Open  the  directory  file  for  input... 

? 

OPEN/READ  INPUT_FILE  'MODEL'. DCT 
» 

LAST^FILE  =  "  " 

I  Loop  through  the  list  of  modified  common  blocks... 

» 

BEGIN: 


READ/END_OF_FILE=END  INPUT_FILE  MODULE_NAME 
MODULE_NAME  -  F$ EDIT ( MODULE_NAME ,  * COMPRE S S , TRIM " ) 
FILE  =  F  $  PARSE { MODULE_NAME ,  ,  r  " NAME  w ) 

IF  FILE  . EQS .  LAST_FILE  THEN  GOTO  BEGIN 

WRITE  SYS  $  OUTPUT  n  "  ,F$EDIT(FILE,  ” LOWERCASE 

OFILE  =  CMNDAT  +  FILE  +  n . "  +  DEXT 

SEARCH/WINDOW=0/OUTPUT= ' OFILE '  'SFILE'  ' FILE ' ”/w LIST 

LAST_FILE  a*  FILE 
GOTO  BEGIN 


END  : 


FINISHED : 

i 

1  Close  and  delete  the  directory  file... 


$  CLOSE  INPUT_FILE 
$  DELETE  ' MODEL '. DCT . * 
$  EXIT 


'  COMMONS ' 
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